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The last decades witnessed a renewal of interest in the Burgers equation. Much activities focused on extensions of the original 
one- dimensional pressureless model introduced in the thirties by the Dutch scientist J.M. Burgers, and more precisely on the 
problem of Burgers turbulence, that is the study of the solutions to the one- or multi-dimensional Burgers equation with random 
initial conditions or random forcing. Such work was frequently motivated by new emerging applications of Burgers model to 
statistical physics, cosmology, and fluid dynamics. Also Burgers turbulence appeared as one of the simplest instances of a 
nonlinear system out of equilibrium. The study of random Lagrangian systems, of stochastic partial differential equations and 
their invariant measures, the theory of dynamical systems, the applications of field theory to the understanding of dissipative 
anomalies and of multiscaling in hydrodynamic turbulence have benefited significantly from progress in Burgers turbulence. 
The aim of this review is to give a unified view of selected work stemming from these rather diverse disciplines. 
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1 Prom interface dynamics to cosmology 

At the end of the thirties, the Dutch scientist J.M. Burg- 
ers [26] introduced a one-dimensional model for pressure- 
less gas dynamics. He was hoping that the use of a simple 
model having much in common with the Navier-Stokes 
equation would significantly contribute to the study of 
fluid turbulence. This model now known as the Burgers 
equation 

dtv + v d x v = vd^v, (1.1) 

has not only the same type of hydrodynamical (or advec- 
tive) quadratic nonlinearity as the Navier-Stokes equa- 
tion that is balanced by a diffusive term, but it also has 
similar invariances and conservation laws (invariance un- 
der translations in space and time, parity invariance, 
conservation of energy and momentum in one dimension 
fori/ = 0). 

Such hopes appeared to be shattered when in the fifties, 
Hopf [67] and Cole [33] showed that the Burgers equation 
can be integrated explicitly. This model thus lacks one 
of the essential properties of Navier-Stokes turbulence: 
sensitivity to small perturbations in the initial data and 
thus the spontaneous arise of randomness by chaotic dy- 
namics. Unable to cope with such a fundamental aspect, 
the Burgers equation then lost its interest in "explain- 
ing" fluid turbulence. 

In spite of this, the Burgers equation reappeared in the 
eighties as the asymptotic form of various nonlinear dis- 
sipative systems. Physicists and astrophysicists then de- 
voted important effort to the understanding of its multi- 
dimensional form and to the study of its random solu- 
tions arising from random initial conditions or a random 



forcing. The goal of this paper is to review selected works 
that exemplify this strong renewal of interest in Burgers 
turbulence. 

The rest of this section is dedicated to the description of 
several physical situations where the Burgers equation 
arises. We will then see in section 2 that in any dimen- 
sion and in the limit of vanishing viscosity the solutions 
to the Burgers equation can be expressed in an explicit 
manner in the decaying case or in an implicit manner in 
the forced case, in terms of a variational principle that 
permits a systematic classification of its various singu- 
larities (shocks and others) and of their local structure 
(normal form). Section 3 is dedicated to the study of the 
decay of the solutions to the one-dimensional unforced 
Burgers equation with random initial data. The multi- 
dimensional decaying problem is discussed in section 4. 
The motivation comes from cosmology where large-scale 
structures can be described in terms of mass transport 
by solutions to the Burgers equation. The basic prin- 
ciples of the forced Burgers turbulence are discussed 
in section 5 where the notions of global minimizer and 
topological shocks are introduced. Section 6 is dedicated 
to the study of the solutions to the periodically kicked 
Burgers equation and their relation with Aubry-Mather 
theory for commensurate- incommensurate phase transi- 
tions. Section 7 reviews various studies of the stochas- 
tically forced Burgers equation in one dimension with a 
particular emphasize on questions that are arising from 
the statistical study of turbulent flows. Finally, section 8 
encompasses concluding remarks and a non-exhaustive 
list of open questions on the problem of the Burgers tur- 
bulence. 

1.1 The Burgers equation in statistical mechanics 

The Burgers equation appears in condensed matter, in 
statistical physics, and also beyond physics in vehicle 
traffic models (see [32], for a review on this topic). When 
a random forcing term is added - usually a white noise 
in time - it is used to describe various problems of in- 
terface deposition and growth (see, for instance, [5]). An 
instance frequently studied is the Kardar-Parisi-Zhang 
(KPZ) model [74]. This continuous version of ballistic 
deposition models accounts for the lateral growth of the 
interface. Let us indeed consider an interface where par- 
ticles deposit with a random flux F that depends both 
on time t and on the horizontal position x. The growth 
of the local height h happens in the direction normal to 
the interface and its time evolution is given by 

dth- i |V/i| 2 = vX7 2 h + F 7 (1.2) 

where the first term of the right-hand side represents the 
relaxation due to a surface tension v. The gradient of 
(1.2) gives the multidimensional Burgers equation 

d t v + v-Vv = vV 2 v-VF, v = -Vh, (1.3) 
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forced by the random potential F. As we will see later, 
shocks generically appear in the solution to the Burgers 
equation in the inviscid limit v — > 0. They correspond 
to discontinuities of the derivative of the height h. The 
KPZ model is hence frequently used to understand the 
appearence of roughness in various interface problems, as 
for instance front propagation in randomly distributed 
forests (see, e.g., [101]). 

The Hopf-Cole transformation Z = exp{h/2v) allows 
rewriting (1.2) as a linear problem with random coeffi- 
cients. 



d t Z = v\t 2 Z 



1 

2v 



(1.4) 



This equation appears in many complex systems, as for 
instance directed polymers in random media [75,22]. 
Indeed the solution Z(x,t) is exactly the partition 
function of an elastic string in the random potential 
(l/2z/) F(x,t), subject to the constraint that its bound- 
ary is fixed at (x, t). Note that here, the time variable t 
is actually a space variable in the main direction of the 
polymer. 



1.2 The adhesion model in cosmology 

The multidimensional Burgers equation has important 
applications in cosmology where it is closely linked to 
what is usually referred to as the Zel'dovich approxi- 
mation [112]. In the limit of vanishing viscosity v — > 
the Burgers equation is known as the adhesion model 
[62]. Right after the decoupling between baryons and 
photons, the primitive Universe is a rarefied medium 
without pressure composed mainly of non-collisional 
dust interacting through Newtonian gravity. The initial 
density of this dark matter fluctuates around a mean 
value p. These fluctuations are responsible for the for- 
mation of the large-scale structures in which both the 
dark non-baryonic matter and the luminous baryonic 
matter concentrate. A hydrodynamical formulation of 
the cosmological problem leads to a description where 
matter evolves with a velocity v, solution of the Euler- 
Poisson equation (see, e.g., [98], for further details). 

In the linear theory of the gravitational instability, that 
is for infinitesimally small initial fluctuations of the den- 
sity field, an instability is obtained with potential dom- 
inant modes (i.e. v — — V 1 ?) and, in the suitable coordi- 
nates, the gravitational interactions can be neglected. In 
1970, Zel'dovich proposed to extend these two proper- 
ties to the nonlinear regimes where density fluctuations 
become important. In this approximation, he also pos- 
tulates that the acceleration is a Lagrangian invariant, 
leading to the formation of caustics. N-body simulations 
however show that the large-scale structures of the Uni- 
verse are much simpler than caustics: they resemble sort 
of thin layers in which the particles tend to be trapped 
(see figure 1(a)). 




(b) 



Fig. 1. (a) Projection of the matter distribution in a slice 
obtained from an N-body simulation by the Virgo consor- 
tium [71]. (b) Composite picture showing the superposition 
of the results of an N-body simulation with the skeleton of 
the results obtained from the adhesion model (from [78] ) . 

It was shown by Gurbatov and Saichev [62] that these 
structures are very well approximated by those obtained 
when constraining the particles not to cross each other 
but to stick together. Even if this mechanism is not col- 
lisional but rather gravitational (probably due to insta- 
bilities at small spatial scales), its effect can be modeled 
by a small viscous diffusive term in the Euler-Poisson 
equation and thus amounts to considering the Burgers 
equation in the limit of vanishing viscosity. 

1.3 A benchmark for hydrodynamical turbulence 

As a nonlinear conservation law, and since its solution 
can be easily known explicitly, the one-dimensional 
Burgers equation frequently serves as a testing ground 
for numerical schemes, and especially for those dedi- 
cated to compressible hydrodynamics. For instance, it 
is a central example for the validation of finite- volumes 
schemes. 

The Burgers equation was also used for testing statisti- 
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cal theories of turbulence. For instance, field theoretical 
methods have frequently been applied to turbulence (see 
[96,102]). These approaches had very little impact until 
recently when they led to significant advances in the un- 
derstanding of intermittency in passive scalar advection 
(see, e.g., [46] for a review). In the past such attempts 
were mostly based on a formal expansion of the non- 
linearity using, for instance, Feynman graphs. Since the 
Burgers equation has the same type of quadratic nonlin- 
earity as the Navier-Stokes equation, such methods are 
applicable in both instances. From this point of view, 
it is important to know answers for Burgers turbulence 
to questions that are generally asked for Navier-Stokes 
turbulence. For instance, Burgers turbulence with a ran- 
dom forcing is the counterpart of the hydrodynamical 
turbulence model where a steady state is maintained by 
an external forcing. The Burgers equation has frequently 
been used as a model where the dissipation of kinetic 
energy remains finite in the limit of vanishing viscosity 
(dissipative anomaly). This allows singling out artifacts 
arising from manipulation that ignore shock waves (see, 
for instance, [51,40]). 

Beyond statistical theory, Burgers turbulence gives a 
simple hydrodynamical training ground for developing 
mathematical tools to study not only Navier-Stokes tur- 
bulence but also various hydrodynamical or Lagrangian 
problems. The forced Burgers equation has recently been 
at the center of studies that allowed unifying different 
branches of mathematics. Mainly used in the past as a 
simple illustration of the notion of entropy (or viscos- 
ity) solution for conservation laws [83,95,85], the Burg- 
ers equation was related in the eighties to the theory 
of Hamiltonian systems developed by Kolmogorov [80] , 
Arnold [2] andMoser [93] (KAM), through the introduc- 
tion of the weak KAM theory [43,47,48]. More recently, 
the study of the solutions to the Burgers equation with a 
random forcing was at the center of a "random" Aubry- 
Mather theory related to random Lagrangian systems 
[38,69]. A particular emphasis on these aspects of Burg- 
ers turbulence is given throughout the present review. 
For the application of the Burgers equation to the prop- 
agation of random nonlinear waves in nondispersive me- 
dia, we refer the reader to the book written by Gurba- 
tov, Malakhov, and Saichev [61]. For a complete state of 
the art on most mathematical apsects of Burgers turbu- 
lence, we refer the reader to the lecture notes by Woy- 
czyhski [110]. 

2 Basic tools 

In this section we introduce various analytical, geomet- 
rical and numerical tools that are useful for constructing 
solutions to the Burgers equation, with and without forc- 
ing, in the limit of vanishing viscosity. All these tools are 
derived from a variational principle that allows writing 
in an implicit way the solution at any time. This varia- 
tional principle leads to a straightforward classification 



of the various singularities that are generically present 
in the solution to the Burgers equation. 

2. 1 Inviscid limit and variational principle 

We consider here the multidimensional viscous Burgers 
equation with forcing 



d t v + (v ■ V)v = vV 2 v - VF(x, t) 



(2.1) 



where x lives on a prescribed configuration space f2 of 
dimension d. For a potential initial condition, v(x, to) = 
— V'I'o(aO) the velocity field remains potential by con- 
struction at any later time, v = — V^, where the poten- 
tial ^ satisfies the equation 



d t y - I |V*| 2 = ^V 2 * + F. 



(2.2) 



Note that if one sets abruptly v = in (2.2), then — 
solves the Hamilton-Jacobi equation associated to the 
Hamiltonian H(q,p) = \p\ 2 + F(q,t). In the unforced 
case, — "3/ is a solution of the Hamilton-Jacobi equation 
associated to the dynamics of free particles. The Hopf- 
Cole transformation [67,33] uses a change of unknown 
^(x, t) = 2v In 6(5;, t). The new unknown scalar field 6 
is solution of the (imaginary-time) Schrodingcr equation 



d t Q = i/V 2 e + ^fe, 

2v 



(2.3) 



with the initial condition 6(5;, to) = exp(^fo(x)/(2v)). 
The solution can be expressed through the Feynman- 
Kac formula 



6(5?, t) = ( cxp 



1 



r*o(^)-r F(W s ,s)ds 



2;/ 



2v 



to 



, (2-4) 



where the brackets (•) denote the ensemble average with 
respect to the realizations of the d-dimensional Brown- 
ian motion W s with variance 2v defined on the configu- 
ration space O and which starts at x at time t. The limit 
v — > is obtained by a classical saddle-point argument. 
The main contribution will come from the trajectories 
W minimizing the argument of the exponential; the ve- 
locity potential can then be expressed as a solution of 
the variational principle 



*(*,*) = -inf [Atf,t ,t) 

7(0 



tt„(7(*o))] , (2-5) 



where the infimum is taken over all trajectories 7 that 
are absolutely continuous (e.g. piece-wise diffcrcntiable) 
with respect to the time variable and that satisfy j(t) = 
x. The action A(j,t n ,t) associated to the trajectory 7 
is defined by 



-4(7, t 



<> = /' 



tf( s )f-Ftf( S ),s) 



ds, (2.6) 
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where the dot stands for time derivative. The kinetic 
energy term |7| 2 /2 comes from the propagator of the d- 

dimensional Brownian motion W . This variational for- 
mulation of the solution to the Burgers equation was 
obtained first by Hopf [67], Lax [83] and Olcinik [95] 
for scalar conservation laws. Its generalization to mul- 
tidimensional Hamilton-Jacobi equations was done by 
Kruzhkov [82] (see also [85]). In the case of a random 
forcing potential F, it was shown by E, Khanin, Mazel 
and Sinai [38] that this formulation is still valid after re- 
placing the action by a stochastic integral. It is also im- 
portant to notice that the variational formulation (2.5) 
in the limit of vanishing viscosity is valid irrespective of 
the configuration space £1 on which the solution is de- 
fined. 

The minimizing trajectories j(-) necessarily satisfy at 
times s < t the Newton (or Eulcr-Lagrange) equation 



2.2 Variational principle for the viscous case 

The derivation of the variational principle (2.5) makes 
use of the Hopf-Cole transformation and of the 
Fcynman-Kac formula. There is in fact another ap- 
proach which also yields a variational formulation of 
the solution to the viscous Hamilton-Jacobi equation 
(2.2). Indeed it turns out that the solution to (2.2) can 
be obtained in the following way. Consider solutions to 
the stochastic differential equation 



(2.11) 



&lu = w(7a, s) ds + \J2v dW ( 



where ft is a stochastic control, that is an arbitrary time- 
dependent velocity field which depends (progressively 
measurably) on the noise W. Limiting ourselves to so- 
lutions satisfying the final condition 73 (t) — x, we can 
write 



7 = -VF(7(s),s), (2.7) 
with the boundary conditions (at the final time t) 



7(i) = x and j(t) = v(x,t). 



(2.8) 



Note that these equations are only valid backward in 
time. Extending them to times larger than t requires 
knowing that the Lagrangian particle will neither cross 
the trajectory of another particle, nor be absorbed by 
a mature shock. This requires global knowledge of the 
solution that satisfies the variational principle (2.5). 



tf (z, t) = - inf (And*, to, t) - (7a(*o))> , (2.12) 

u 

where the brackets (•) now denote average with respect 
to W s and the action is given by 



Afi("fji, to, t) 



1 



u(s)\ 2 -F(Us),s) 



ds. (2.13) 



It is obvious that this variational principle gives (2.5) in 
the inviscid limit v — ► 0. Note that this approach has the 
advantage to be applicable not only to Burgers dynamics 
but to any convex Lagrangian (see [50,58]). 



When the forcing term is absent from (2.1), it is easily 
checked that the variational principle reduces to 



max * (^o) - 

x 



\x — xq\ 2 
2t 



(2.9) 



where the maximum is taken over all initial positions 
xq in the configuration space O. The Euler-Lagrangc 
equation takes then the particularly simple form 



7 = 0, i.e. x = x Q + tv (x ), 



(2.10) 



which simply means that the initial velocity is conserved 
along characteristics. 

Typically there exist Eulerian locations x where the min- 
imum in (2.5) - or the maximum in (2.9) in the unforced 
case - is reached for several different trajectories 7- Such 
locations correspond to singularities in the solution to 
the Burgers equation. After their appearance, the veloc- 
ity potential contains angular points corresponding to 
discontinuities of the velocity field v = —V*. 



2.3 Singularities of Burgers turbulence 

The singularities appearing in the course of time play 
an essential role in understanding various aspects of the 
statistical properties in the inviscid limit. The shocks - 
discontinuities of the velocity field - and other singular- 
ities, such as preshocks, generally not associated to dis- 
continuities, are often responsible for non-trivial univer- 
sal behaviors. In order to understand the contribution 
of each kind of singularities, it is first important to know 
in a detailed manner their genericity and their type. 

As we have seen in the previous section, the potential 
solutions to the multidimensional Burgers equation can 
be expressed in the inviscid limit in terms of the vari- 
ational principle (2.5) (that reduces to (2.9) in the un- 
forced case). There typically exist Eulerian locations x 
where the minimum is cither degenerate or attained for 
several trajectories. A co-dimension can be associated 
to such points by counting the number of relations that 
are necessary to determine them. The singular locations 
of co-dimension c form manifolds of the Eulerian space- 
time with dimension (d — c) . The singularities with the 
lower co-dimension are the shocks corresponding to the 
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Eulcrian positions where two different trajectories min- 
imize (2.5); they form Eulerian manifolds of dimension 
(d— 1): in one dimension the shocks are isolated points, 
in two dimensions they are lines, in three dimensions 
surfaces, etc. There also exist Eulerian manifolds with 
three different minimizing trajectories. In one dimen- 
sion, they are isolated space-time events corresponding 
to the merger of two shocks. In two dimensions, they are 
triple points where three shock lines meet. In three di- 
mensions they are filaments corresponding to the inter- 
section of three shock surfaces. There also exist Eulerian 
locations where the minimum in (2.5) is reached for four 
different trajectories, etc. 
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Fig. 2. Typical aspect of the singularities present at a fixed 
time in the solution for (a) d = 2 and (b) d — 3. 

The generic form of such singularities and their typi- 
cal metamorphoses occurring in the course of time were 
studied in details and classified for d = 2 and d — 3 by 
Arnold, Baryshnikov and Bogaevsky in the Appendix 
of [62] and in a more detailed paper by Bogaevsky [17]. 
This classification is based on two criteria: (i) the num- 
ber of trajectories minimizing (2.5) and (ii) the multi- 
plicity of each of these minima. The shocks correspond- 
ing to locations with two distinct minimizers are hence 
denoted by A\. At a fixed time, the A\ singularities are 
discrete points in one dimension. In two dimensions (see 
figure 2(a)) they form curve segments with extremities 
that can be either triple points A\ or isolated termina- 



tion points of the type A3 corresponding to a degenerate 
minimum. In three dimensions (see figure 2(b)) the sin- 
gular manifold is formed by shock surfaces of A\ points. 
The boundaries of these surfaces are either made of de- 
generate A 3 points or of triple lines made of A\ points. 
The triple lines intersect at isolated A\ points or inter- 
sect shock boundaries at particular singularities called 
AiA 3 where the minimum is attained in two points, one 
of which is degenerate. 

It is important to remark here that degenerate singu- 
larities (of the type A 3 or of higher orders A 5 , At, etc.) 
introduce in the solution points where the velocity gra- 
dients becomes arbitrarily large. This is not the case of 
the A™ singularities which correspond to discontinuities 
of the velocity but are associated to bounded values of 
its gradients. As we will see in sections 4 and 7, these 
degenerate singularities are responsible for an algebraic 
behavior of the probability density function of velocity 
gradients, velocity increments and of the mass density. 




Fig. 3. Illustration of the similarities between the singular 
manifold in space time for d — 1 and at fixed time for d = 2 
(b). The two manifolds contain the same type of singularities 
with the same co-dimensions. The restrictions on the possible 
metamorphoses in dimension d = 1 are the following: a point 
of the type A3 can only exist at the bottom extremity of 
a shock trajectory; the A? points necessarily correspond to 
the merger of two shocks; shock trajectories cannot have a 
horizontal tangent. 

The singularities with co-dimensions (d+ 1) gcnerically 
appear in the solution at isolated times. They corre- 
spond to instantaneous changes in the topological struc- 
ture of the singular manifold, called metamorphoses and 
can be also classified (see [17]). In one dimension, there 
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are two generic metamorphoses: shock formations (the 
preshocks) corresponding to a specific space-time loca- 
tion where the minimum is degenerate (A 3 singularities) 
and shock mergers associated to space-time positions 
where the minimum is attained for three different tra- 
jectories (A\ singularities). We see that some of the sin- 
gularities generically present in two dimensions appear 
at isolated times in three dimensions. Actually, all the 
singularities generically present in dimension (d+ 1) ap- 
pear in dimension d on a discrete set of space time, that 
is at isolated positions and instants of time. However, it 
has been shown in [17] that the irreversible dynamics of 
the Burgers equation restricts the set of possible meta- 
morphoses. The admissible metamorphoses are charac- 
terized by the following property: after the bifurcation, 
the singular manifold must remain locally contractiblc 
(homotopic to a point in the neighborhood of the Eu- 
lerian location of the metamorphosis). This topological 
restriction is illustrated for the one-dimensional case in 
figure 3. Note that this constraint actually holds for all 
solutions to the Hamilton-Jacobi equation in the limit 
of vanishing viscosity, as long as the Hamiltonian is a 
convex function. 

In order to determine precisely how all these singularities 
contribute to the statistical properties of the solution, it 
is important to know the local structure of the velocity 
(or potential) field in their vicinity. Various normal forms 
can be obtained from the multiplicity of the minimum 
in the variational formulation of the solution (2.5). In 
the case without forcing, they can be obtained from a 
Taylor expansion of the initial velocity potential. This 
will be used in next section to determine the tail of the 
probability distribution of a mass density field advected 
by a velocity solution to the Burgers equation. 



2.4 Remarks on numerical methods 

All the traditional methods used to solve equations of 
fluid dynamics, or more generally any partial differen- 
tial equations, can be used to obtain the solutions to the 
Burgers equation. However, as we have seen above, the 
solution typically has singularities (discontinuities of the 
velocity) in the limit of vanishing viscosity. Hence meth- 
ods which rely on the smoothness of the solution require 
a non-vanishing viscosity, which is introduced either in 
an explicit way to ensure stability (as, e.g., for pseudo- 
spectral methods) or in an implicit way through the 
discretization procedure (as for finite-differences meth- 
ods). In both cases the value of the viscosity is deter- 
mined from the mesh size and, even in one dimension, 
their uses might be very disadvantageous. We will now 
demonstrate various numerical methods that allow ap- 
proximating the solutions to the Burgers equation di- 
rectly in the limit of vanishing viscosity v — > 0. 



2.4-1 Finite volumes 

The one-dimensional Burgers equation with no forcing is 
a scalar conservation law. Its entropic solutions (or vis- 
cosity solutions) can thus be approximated numerically 
by finite- volume methods. Instead of constructing a dis- 
crete approximation of the solution on a grid, such meth- 
ods consist in considering an approximation of its mean 
value on a discrete partitioning of the system into finite 
volumes. One then needs to evaluate for each of these 
volumes the fluxes exchanged with each of its neighbors. 
Various approximations of these fluxes were introduced 
by Godunov, Roe, and Lax and Wendroff (see, e.g., [35], 
Vol. 3, for a review). These methods require to dicretize 
both space and time. The time step being then related 
to the spatial mesh size by a Courant-Friedrichs-Lewy 
type condition. Thus to integrate the equation during 
times comparable to one eddy turnover time, they re- 
quire a computational time 0(N 2 ) where N is the res- 
olution. As we now show there actually exist numerical 
schemes that allow constructing the solution to the de- 
caying Burgers equation for arbitrary times without any 
need to compute the solution at intermediate times. 

2-4-2 Fast Legendre transform 

As we have seen in section 2.1, the solution to the un- 
forced Burgers equation is given explicitly by the varia- 
tional principle (2.9). A method based on the idea of us- 
ing this formulation together with a monotonicity prop- 
erty of the Lagrangian map x — > x — X(xo,t) was 
given in [94]. It is called the fast Legendre transform 
whose principles were already sketched in [23] . Both Eu- 
lerian and Lagrangian positions are discretized on reg- 
ular grids. Then, for a fixed Eulerian location x« on 
the grid, one has to find the corresponding Lagrangian 

coordinate Xq^ maximizing (2.9). A naive implementa- 
tion would require 0(N£ N£) operations if the Eulerian 
and the Lagrangian grids contain and points re- 
spectively. Actually the number of operations can be re- 
duced to 0((Ne lnAL) d ) by using the monotonicity of 
the Lagrangian map, that is the fact that for any pair 
of Lagrangian positions x^ and Xq \ one has at any 

time [X(xf\t) - X(xf\t)] ■ (xf ] - x^) > 0. In the 
case of orthogonal grids, this property allows perform- 
ing the maximization by exploring along a binary tree 
the various possibilities; thus the number of operations 
is reduced to In iV L for each of the N& positions on the 
Eulerian grid. Such algorithms give access to the solu- 
tion not only directly in the limit of vanishing viscosity 
but also by jumping directly from the initial time to an 
arbitrary time. 

This method can also be used for the forced Burgers 
equation, approximating the forcing by a sum of im- 
pulses at discrete times and letting the solution decay 
between two such kicks. This gives an efficient algorithm 
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for the forced Burgers equation directly applicable in the 
limit of vanishing viscosity. 

2.4-3 Particle tracking methods 

In one dimension, Lagrangian methods can be imple- 
mented in a straightforward manner after noticing that 
particles cannot cross each other and that it is advis- 
able to track not only fluid particles but also shocks (see, 
e.g., [6]). Lagrangian methods can in principle be used 
to solve the Burgers equation in any dimension. How- 
ever the shock dynamics is meaningful only for poten- 
tial solutions. Outside the potential framework almost 
nothing is known about the construction of the solution 
beyond the first crossing of trajectories. In the potential 
case, a particle method can be formulated by choosing 
to represent the solution in the position-potential (x, 'J) 
space instead of the position-velocity (x, v) space. An 
idea in two dimensions, which was not yet implemented, 
consists in considering a meshing of the hyper-surface 
defined by the velocity potential. If such a mesh contains 
only triple points, such singularities are preserved by the 
dynamics and can be tracked using the results discussed 
below in section 4.2 and by checking at all time steps in 
an exhaustive manner at all the metamorphoses encoun- 
tered by triple points. 



3. 1 Geometrical constructions of the solution 
3.1.1 The potential Lagrangian manifold 

The variational formulation of the solution (3.2) has 
a simple geometrical interpretation in the position- 
potential space (x, of dimension d + 1. Indeed, con- 
sider the ci-dimensional manifold parameterized by the 
Lagrangian coordinate xo and defined by 

!x = x - t V* (£o) 
t n (3-3) 

* = * (x )--|V* (f )| 2 . 

The first line corresponds to the position where the gra- 
dient of the argument of the maximum function in (3.2) 
vanishes while the second line is just its argument eval- 
uated at the maximum. For a sufficiently regular initial 
potential \Po (at least twice differentiable) and for suf- 
ficiently small times, equation (3.3) unambiguously de- 
fines a single- valued function ^(x, t). However, there ex- 
ists generically a time at which the manifold is folding. 
Figure 4(a) (upper) shows in one space dimension the 
typical shape of the Lagrangian manifold defined by (3.3) 
after the critical time t*. For some Eulerian positions x, 
there is more than one branch and cusps are present at 
Eulerian locations where the number of branches change. 



3 Decaying Burgers turbulence 

We focus in this section on the solutions to the d- 
dimensional unforced potential Burgers equation 

d t v + v-Vv = vV 2 v, v{x,Q)=v {x) = -V^o{x). (3.1) 

As showed in section 2.1, the solution can be expressed 
in the limit of vanishing viscosity v —* in terms of a 
variational principle that relates the velocity potential 
at time t to its initial value: 

*(f , t) = max U {x Q ) - |f ~f°H . (3.2) 

x a \ At J 

The next subsection describes several geometrical con- 
structions of the solution that are helpful to determine 
various statistical properties of the decaying prob- 
lem (3.1). This is illustrated in subsections 3.2 and 3.3 
which are devoted to the study of the decay of smooth 
homogeneous and of Brownian initial data, respectively. 

The study of the solutions to the Burgers equation trans- 
porting a density field is of particular interest in the ap- 
plication of the Burgers equation in cosmology within 
the framework of the adhesion model. This question will 
be discussed in section 4. 



The situation is very similar in higher dimensions as il- 
lustrated for d = 2 in figure 4(b). Clearly from the varia- 
tional principle (2.9), the correct solution to the inviscid 
Burgers equation is obtained by taking the maximum, 
that is the highest branch. The velocity potential is by 
construction always continuous but it contains angular 
points corresponding to discontinuities of the velocity 
v = — V*. Such singularities are located at Eulerian lo- 
cations where the maximum in (2.9) is degenerate and 
attained for different x . As already discussed in sec- 
tion 2.3 the different singularities appearing in the solu- 
tions can be classified in any dimension. 

Below we describe other geometrical constructions of the 
solutions to the decaying Burgers equation in the limit 
of vanishing viscosity that are based on the variational 
principle (2.9). 

3.1.2 The velocity Lagrangian manifold 

In one dimension, when the velocity field is always po- 
tential, the method based on the study of the poten- 
tial manifold in the (x, space described above can 
be straightforwardly extended to the position-velocity 
phase space. Consider the Lagrangian manifold defined 

by 

x = x a -tv (x ) ^ ^ 

V = Vq(xq). 
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Maxwell 




(b) 

Fig. 4. (a) Lagrangian manifold for d = 1 in the (a;, \t) 
plane (upper) and in the (x, v) plane (lower); the heavy lines 
correspond to the correct Eulerian solutions, (b) Lagrangian 
manifold in the (x, \t) space for d = 2. 

The regular parts of the graph of the solution are nec- 
essarily contained in this manifold. However, for times 
larger than i*, folding appears and the naive solution 
would be multi-valued. To construct the true solution 
one should find a way to choose among the different 
branches. In one dimension, there is a simple relation 
between the potential Lagrangian manifold in the (x, \I/) 
plane and those of the (x,v) plane defined by (3.4): 
the potential manifold is obtained by taking the "multi- 
valued integral" that can be defined by transforming the 
spatial integral into an integral with respect to the arc 
length. The maximum representation (2.9) implies that 
the velocity potential is continuous. Hence a shock cor- 
responds to an Eulerian position x where two points be- 
longing to different branches define equal areas in the 
(x, v) plane. In the case of a single loop of the manifold, 
this is equivalent to applying the Maxwell rule to deter- 
mine the shock position (see figure 4(a) - lower). This 
construction of the solution can become rather involved 
as soon as the number of shocks becomes large or that 
several mergers have taken place. For the moment there 



shock x 
interval 

(a) 



Ruled surface 




(b) 

Fig. 5. Convex hull construction in terms of the Lagrangian 
potential (a) for d = 1 and (b) for d — 2. 

is no generalization to dimension higher than one of this 
Maxwell rule construction of the solution. For such an 
extension, one needs to develop a geometrical framework 
to describe the Lagrangian manifold in the (x, v) space. 
Such approaches would certainly shed some light on the 
problem of constructing non-potential solutions to the 
Burgers equation in the limit of vanishing viscosity. 

3.1.3 The convex hull of the Lagrangian potential 

Another geometrical construction of the solution, which 
is valid in any dimension makes use of the Lagrangian 
potential 

I ~* |2 

$(*o,t)=t*o(*o)-^-. (3.5) 

Clearly, the negative gradient of the Lagrangian poten- 
tial gives the naive Lagrangian map 

X(x , t) = -Vx ${x*o,t) =x + tv (x ), (3.6) 

that is satisfied by Lagrangian trajectories as long as 
they do not enter shocks. The maximum formulation of 
the solution (2.9) can be rewritten as 

|f|2 

W(x, t) + LJ- = max($(.f ) + x ■ x), (3.7) 

2 x 
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which represents the potential as, basically, a Legcndre 
transform of the Lagrangian potential. An important 
property of the Legendre transform is that the right- 
hand side, of (3.7) is unchanged if the Lagrangian po- 
tential $ is replaced by its convex hull, that is the in- 
tersection of all the half planes containing its graph. In 
other terms, the convex hull $ c of the Lagrangian po- 
tential $ is defined as <& c (x ,i) = infg(x ), where the 
infimum is taken over all convex functions g satisfying 
g(-) > This is illustrated in one dimension in fig- 

ure 5(a) which shows both regular points (Lagrangian 
points which have not fallen into a shock) and one shock 
interval, situated below the segment which is a part of 
the convex hull. In two dimensions, as illustrated in fig- 
ure 5(b), the convex hull is typically formed by regular 
points, by ruled surfaces, and by triangles which corre- 
spond, to the regular part of the velocity field, the shock 
lines, and the shock nodes, respectively. 

Note that in one dimension, there exists an equivalent 
construction which is directly based on the Lagrangian 
map xo i— > X(xo,t) defined by (3.6). Working with the 
convex hull is equivalent to the Maxwell rule applied 
to the non-invertible regions of the Lagrangian map. A 
shock corresponds to a whole Lagrangian interval having 
a single point as an Eulerian image. One then talks about 
a Lagrangian shock interval. 

3.1.4 The paraboloid construction 

Finally, the maximum representation (3.7) leads in a 
straightforward way to another geometrical construction 
of the solution. As illustrated in figure 6 in both one and 
two dimensions, a paraboloid with apex at x and ra- 
dius of curvature proportional to t is moved down in the 
(xq, ^0) space until it touches the surface defined by the 
initial velocity potential ^0 at the Lagrangian location 
associated to x. The location xq where the paraboloid 
touches the graph of the potential is exactly the pre- 
image of x. If it touches simultaneously at several loca- 
tions, a shock is located at the Eulerian position x. One 
constructs in this way the inverse Lagrangian map. 

3.2 Kida's law for energy decay 

An important issue in turbulence is that of the law of 
decay at long times when the viscosity is very small. 
Before turning to the Burgers equation it is useful to re- 
call some of the features of decay for the incompressible 
Navier-Stokes case. It is generally believed that high- 
Reynolds number turbulence has universal and non- 
trivial small-scale properties. In contrast, large scales, 
important for practical applications such as transport of 
heat or pollutants, are believed to be non-universal. This 
is however so only for the toy model of turbulence main- 
tained by prescribed large-scale random forces. Very 
high-Reynolds number turbulence, decaying away from 
its production source, and far from boundaries can relax 




(b) 

Fig. 6. Paraboloid construction of solution for (a) d = 1 and 
(b) d = 2. 

under its internal nonlinear dynamics to a (self-similarly 
evolving) state with universal and non-trivial statistical 
properties at all scales. Von Karman and Howarth [109], 
investigating the decay for the case of high-Reynolds 
number homogeneous isotropic three-dimensional tur- 
bulence, proposed a self-preservation (self-similarity) 
ansatz for the spatial correlation function of the ve- 
locity: the functional shape of the correlation function 
remains fixed, while the integral scale L(t) grows in 
time and the mean kinetic energy E(t) = V 2 (t) decays, 
both following power laws; there are two exponents 
which can be related by the condition that the energy 
dissipation per unit mass |-E(t)| should be proportional 
to V 3 /L. But an additional relation is needed to actu- 
ally determine the exponents. The invariance in time of 
the energy spectrum at low wavenumbers, known as the 
"permanence of large eddies" [53,84,63] can be used to 
derive the law of self-similar decay when the initial spec- 
trum Eo(k) cx k n at small wavenumbers k with n below 
a critical value equal to 3 or 4, the actual value being 
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Fig. 7. Snapshot of solution of decaying Burgers turbulence 
at long times when spatial periodicity is assumed. 

disputed because of the "Gurbatov phenomenon" (see 
the end of this section). One then obtains a law of decay 
E(t) oc t-2(n+i)/(3+n)_ (Kolmogorov [79] proposed a law 
of energy decay V 2 (t) oc i~ 10 / 7 , which corresponds to 
n = 4 and used in its derivation the so-called "Loitsyan- 
sky invariant", a quantity actually not conserved, as 
shown by Proudman and Reid [100].) When the initial 
energy spectrum at low wavenumbers goes to zero too 
quickly, the permanence of large eddies cannot be used, 
because the energy gets backscattered to low wavenum- 
bers by nonlinear interactions. For Navier-Stokes tur- 
bulence the true law of decay is then known only within 
the framework of closure theories (see, e.g., [84]). 

For one-dimensional Burgers turbulence, many of the 
above issues are completely settled. First, we observe 
that the problem of decay is quite simple if spatial peri- 
odicity is assumed. Indeed, all the shocks appearing in 
the solution will eventually merge into a single shock per 
period, as shown in figure 7. The position of this shock 
is random and the two ramps have slope 1 /t, as is easily 
shown using the parabola construction of subsection 3.1. 
Hence, the law of decay is simply E(t) oc t~ 2 . Nontrivial 
laws of decay are obtained if the Burgers turbulence 
is homogeneous in an unbounded domain and has the 
"mixing" property (which means, roughly, that correla- 
tions are vanishing when the separation goes to infinity). 
The number of shocks is then typically infinite but their 
density per unit length decreases in time because shocks 
are constantly merging. The E{t) oc i- 2 ("+ 1 )/( 3 +™) l aw 
mentioned above can be derived for Burgers turbulence 
from the permanence of large eddies when n < 1 [63]. 
For n = 0, this i~ 2 / 3 law was actually derived by Burg- 
ers himself [27]. 

The hardest problem is again when permanence of 
large eddies does not determine the outcome, namely 
for n > 1. This problem was solved by Kida [77] (see 
also [51,61,63]). 

We now give some key ideas regarding the derivation of 
Kida's law of energy decay. We assume Gaussian, homo- 
geneous smooth initial conditions, such that the poten- 
tial is homogeneous. Note that a homogeneous function 
is not, in general, the derivative of another homogeneous 
function. Here this is guaranteed by assuming that the 



Fig. 8. An initial potential which is everywhere below the 
parabola x /(2t) + The probability of such events gives 
the cumulative probability to have a potential at time t less 
than *. 

initial spectrum of the kinetic energy is of the form 

E (k) oc k n for k -» with n > 1 . (3.8) 

This condition implies that the mean square initial po- 
tential J k~ 2 E (k) dk has no infrared ( small- fc) diver- 
gence (the absence of an ultraviolet divergence is guar- 
anteed by the assumed smoothness). 

A very useful property of decaying Burgers turbulence, 
which has no counterpart for Navier-Stokes turbulence, 
is the relation 

E(t) = | <*> , (3-9) 

which follows by taking the mean of the Hamilton 
Jacobi equation for the potential (2.2) in the absence of 
viscosity and of a driving force. Hence, the law of energy 
decay can be obtained from the law for the mean po- 
tential. The latter can be derived from the cumulative 
probability of the potential which, by homogeneity, does 
not depend on the position. By (2.9), its expression at 
x = is 

Prob{Pot.<*} = Prob|va;o, *o(zo) < t| + *}- ( 3 - 10 ) 

Expressed in words, we want to find the probability 
that the initial potential does not cross the parabola 
£o/(2t) + ^ (see figure 8). Since, at large times t, the 
relevant \& is going to be large, the problem becomes 
that of not crossing a parabola with small curvature and 
very high apex. The crossings, more precisely the up- 
crossings, are spatially quite rare and, for large t, form 
a Poisson process [92] for which 

Prob. no crossing ~ e~^ N ^\ (3-11) 

where (N(t)} is the mean number of up-crossings. By 
the Rice formula (a consequence of the identity S(Xx) = 

<w(t » -(£>'<•"•<«>>-•> S H (S))' <3 - 12 > 
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where H is the Hcavisidc function and 



m(x ) = ^o{x ) - 



(3.13) 



Since ^o^o) is Gaussian, the right-hand side of (3.12) 
can be easily expressed in terms of integrals over the 
probability densities of "^oO^o) and of d^oOcoVdxo (as 
a consequence of homogeneity these variables are uncor- 
rected and, hence, independent). The resulting integral 
can then be expanded by Laplace's method for large t, 
yielding 



(AT(i))~t 1/2 *- 1/2 e-* 2 , t^oo. 



(3.14) 



When this expression is used in (3.11) and the result is 
differentiated with respect to \I/ to obtain the probability 
density function (PDF) of p(^>), the latter is found to 
be concentrated around = (hit) 1 / 2 . It then follows 
that, at large times, we obtain Kida's log-corrected 1/t 
law for the energy decay 



(*) ~ (lnt) 1/2 , E(t) 



L{t) 



t 
hit 



1/4 

. (3.15) 



The Eulerian solution, at large times, has the ramp 




Fig. 9. The Eulerian solution at large times t. The ramps 
have slope 1/t. In time-independent scales, the figure would 
be stretched horizontally and squeezed vertically by a factor 
proportional to t. 

structure shown in hgure 9 with shocks of typical 
strength V(t) = E x / 2 {t), separated by a distance L(t). 

The fact that Kida's law is valid for any n > 1, and not 
just for n > 2 as thought originally, gives rise to an inter- 
esting phenomenon now known as the "Gurbatov effect" : 
if 1 < n < 2 the large-time evolution of the energy spec- 
trum cannot be globally self-similar. Indeed, the perma- 
nence of large eddies, which is valid for any n < 2 dic- 
tates that the spectrum should preserve exactly its initial 
C n k n behavior at small wavenumbers k, with a constant- 
in-time C n . Global self-similarity would then imply a 
£-2(n+i)/(3+n) j aw f or £ nc ener gy decay, which would 

contradict Kida's law. Actually, as shown in [63], there 
are two characteristic wavenumbers with different time 
dependences, the integral wavenumber fei,(t) ~ (i(t)) -1 
and a switching wavenumber k s {t) below which holds 
the permanence of large eddies. It was shown that the 
same phenomenon is present also in the decay of a pas- 
sive scalar [45] . Whether or not a similar phenomenon is 



present in three-dimensional Navier-Stokes incompress- 
ible turbulence, or even in closure models, is a contro- 
versial matter [44,97]. 

For decaying Burgers turbulence, if we leave aside the 
Gurbatov phenomenon which does not affect energy- 
carrying scales, the following may be shown. If we rescale 
distances by a factor L(t) and velocity amplitudes by a 
factor V(t) — E 1 / 2 ^) and then let t — > oo, the spatial 
(single-time) statistical properties of the whole random 
velocity field become time-independent. In other words, 
there is a self-similar evolution at large times. Hence, 
dimcnsionlcss ratios such as the velocity flatness 



F(t) 



i(v 2 (t))} 2 



(3.16) 



have a finite limit as t — > oo. A similar property holds 
for the decay of passive scalars [28]. We do not know if 
this property holds also for Navier-Stokes incompress- 
ible turbulence or if, say, the velocity flatness grows with- 
out bounds at large times. 

3.3 Brownian initial velocities 

Initial conditions in the Burgers equation that are Gaus- 
sian with a power-law spectrum oc k~ a have been fre- 
quently studied because they belong in cosmology to the 
class of scale-free initial conditions (see [98,34]). We con- 
sider here the one-dimensional case with Brownian mo- 
tion as initial velocity, corresponding to a = 2. 

Brownian motion is continuous but not diffcrcntiablc; 
thus, shocks appear after arbitrarily short times and are 
actually dense (see figure 10(a)). Numerically supported 
conjectures made in [104] have led to a proof by Sinai 
[105] of the following result: in Lagrangian coordinates, 
the regular points, that is fluid particles which have not 
yet fallen into shocks, form a fractal set of Hausdorff 
dimension 1/2. This implies that the Lagrangian map 
forms a Devil's staircase of dimension 1/2 (see figure 11). 
Note that when the initial velocity is Brownian, the La- 
grangian potential has a second space derivative delta- 
correlated in space; this can be approximately pictured 
as a situation where the Lagrangian potential has very 
wild oscillations in curvature. Hence, it is not surprising 
that very few points of its graph can belong to its convex 
hull (see figure 10(b)). 

We will now highlight some aspects of Sinai's proof of 
this result. The idea is to use the construction of the 
solution in terms of the convex hull of the Lagrangian 
potential (see section 3.1), so that regular points are ex- 
actly points where the graph of the Lagrangian poten- 
tial coincides with its convex hull. For this problem, the 
Hausdorff dimension of the regular points is also equal 
to its box-counting dimension, which is easier to deter- 
mine. One obtains it by finding the probability that a 
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(a) 




(b) 



Fig. 10. Snapshot of the solution resulting from Brownian 
initial data in one dimension, (a) Velocity profile at initial 
time t = and at time t = 1; notice the dense proliferation 
of shocks, (b) Lagrangian potential together with its convex 
hull. 



hull. In other words, one looks for points, such as R, with 
the property that the graph of $ lies below its tangent 
at R (sec figure 12). Following Sinai, this can be equiva- 
lently formulated by the box construction with the fol- 
lowing constraints on the graph: 
Left: graph of the potential below the half line T_ , 
Right: graph of the potential below the half line T + , 



1 



Box: < 



enter (AF) with a slope larger 
than that ofT_ by 0(£^ 2 ) 

exit (CB) with a slope less than 
that of T+ by Oil 1 ' 2 ) 

cross (FC) and stay below (ED). 
It is obvious that such conditions ensure the existence of 
at least one regular point, as seen by moving (ED) down 
parallel to itself until it touches the graph. Note that A 
and the slope of (AB) are prescribed. Hence, one is cal- 
culating conditional probabilities; but it may be shown 
that the conditioning is not affecting the scaling depen- 
dence on £. 




Fig. 12. The box construction used to find a regular point R 
within a Lagrangian interval of length I (from [105,107]). 




As the Brownian motion v§(x§) is a Markov process, the 
constraints Left, Box and Right are independent and 
hence, 

preg. ^ ^ Probjregular point in interval of length £} 
= Prob{Lcft} x Prob{Box} x Prob{ Right} . (3.17) 

The sizes of the box were chosen so that Prob {Box} is 
independent of £: 



Prob {Box} 



(3.18) 



Fig. 11. The Lagrangian map looks like a Devil's staircase: 
it is constant almost everywhere, except on a fractal Can- 
tor-like set (from [107]). 



Indeed, Brownian motion and its integral have scaling 
exponent 1/2 and 3/2, respectively, and the problem 
with £ <C 1 can be rescaled into that with £ = 1 without 
changing probabilities. 



small Lagrangian interval of length £ contains at least 
one regular point which belongs simultaneously to the 
graph of the Lagrangian potential $ and to its convex 



It is clear by symmetry that Prob{Left} and Prob{Kight} 
have the same scaling in I. Let us concentrate on 
Prob{ Bight}. We can write the equation for the half 
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line r + in the form 

r + :.T ^$(x^) + « 3 / 2 
d$ 



dx 



K')+7^ 1/2 



(3.19) 



where 7 and S are positive O(l) quantities. Hence, intro- 
ducing a = xq — x'q, the condition Right can be written 
to the leading order as 



/' 

Jo 



v (x ) + j£ 1/2 ] dx + 8£ 3 ' 2 + y > 0, (3.20) 



for all a > 0. By the change of variable a — [3£ and use of 
the fact that the Brownian motion has scaling exponent 
1/2, one can write the condition Right as 



Jo 



{v (x ) + 7) dx > -6, for all a e [0, fr 1 ]. (3.21) 



Without affecting the leading order, one can replace the 
Brownian motion by a stepwise constant random walk 
with jumps of ±1 at integer xo's. The integral in (3.21) 
has a simple geometric interpretation, as highlighted in 
figure 13, which shows a random walk starting from 
the ordinate 7 and the arches determined by succes- 
sive zero-passings. The areas of these arches are denoted 
<S*j Si, ...S n , 



Wd+Y 




S^. S[ St S ^ S 5 S 



Fig. 13. The arches construction which uses the zero-passings 
of a random walk to estimate the integral of Brownian motion 
(from [105,107]). 



It is easily seen that 

Prob{Kight} - Prob{S*i > 0, Si + S 2 > 0, . . . 

Si + • • • + S n > }, (3.22) 

where n = 0(£~ x / 2 ) is the number of zero-passings of 
the random walk in the interval [0, The probability 
(3.22) can be evaluated by random walk methods (see, 
e.g. ,[49], Chap. 12, section 7), yielding 



Prob{ Right} ~ Prob{n first sums > 0} 



(3.23) 



By (3.17), (3.18) and (3.23), the probability to have a 
regular point in a small interval of length £ behaves as 



I 1 / 2 when I — ► 0. Thus, the regular points have a box- 
counting dimension 1/2. 

This rigorous result on the fractal dimension of regular 
points served as a basis in [4] for a proof of the bifrac- 
tality of the inverse Lagrangian map when the initial 
velocity is Brownian. Namely, the moments M q (£) = 
{(xo(x + £) — xo(x))) behave as £ Tq at small separation 
£ and the exponents T q experience the phase transition 



T q = 2q for q < 1/2 
T q = 1 for q > 1/2 



(3.24) 
(3.25) 



At the moment, this is the only rigorous result on the 
bifractal nature of the solutions to the Burgers equation 
in the case of non-differentiable initial velocity. In par- 
ticular, the case of fractional Brownian motion is still 
opened. 

4 Transport of mass in the Burgers/adhesion 
model 

In the cosmological application of the Burgers equation, 
i.e. for the adhesion model, it is of particular interest 
to analyze the behavior of the density of matter, since 
the large-scale structures are characterized as regions 
where mass is concentrated. This is done by associat- 
ing to the velocity field v solution to the d-dimensional 
decaying Burgers equation (3.1), a continuity equation 
for the transport of a mass density field p. In Eulerian 
coordinates, the mass density p satisfies 



dtp + V • (pv) = , p(x,0) = p {x). 



(4.1) 



A straightforward consequence of (4.1) and of the for- 
mulation of Burgers dynamics in terms of characteristics 
X(xo, t) is that, at the Eulerian locations where the La- 
grangian map is invertible, the mass density field p can 
be expressed as 



p(x, t) = 



Po{xo) 



where X (xq ,t)=x, 



J(xo,t) 

and J(x ,t)^det\(dX l )/(dx 1 ) . 



(4.2) 



Large but finite values of the density will be reached 
at locations where the Jacobian J of the Lagrangian 
map becomes very small. As we will see in section 4.1, 
they contribute a power-law behavior in the tail of the 
probability density function of p. 

The expression (4.2) is no more valid when the Jacobian 
vanishes (inside shocks) . Then the density field becomes 
infinite and mass accumulates on the shock. We will see 
in section 4.2 that the evolution of the mass inside the 
singularities of the solution can be obtained as the v — > 
limit of the well-posed viscous problem. Finally, we 
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will discuss in section 4.3 sonic of the applications of 
the Burgers equation to cosmology, and in particular 
how, assuming the dynamics of the adhesion model, the 
question of reconstruction of the early Universe from 
its present state can be interpreted as a convex optimal 
mass transportation problem. 

4-1 Mass density and singularities 

We give here the proof reported in [54] that in any dimen- 
sion large densities are localized near "kurtoparabolic" 
singularities residing on space-time manifolds of co- 
dimension two. In any dimension, such singularities con- 
tribute universal power-law tails with exponent —7/2 
to the mass density probability density function (PDF) 
p{p), provided that the initial conditions are smooth. 

In one dimension, the mass density at regular points can 
be written as 



where 



p(X(x ,t),t) = 



PoOo) 



l-t[(d^ )/(dxt)} 



(4.3) 



We suppose here that the initial density p is strictly 
positive and that both p and ^ are sufficiently regu- 
lar statistically homogeneous random fields. Large val- 
ues of p{x, t) are obtained in the neighborhood of La- 
grangian positions with a vanishing Jacobian, i.e. where 
A 2 '^q{xq)/Ax^ = 1/t. Once mature shocks have formed, 
the Lagrangian points with vanishing Jacobian are in- 
side shock intervals and thus not regular. The only points 
with a vanishing Jacobian that are at the boundary of 
the regular points are obtained at the preshocks, that is 
when a new shock is just born at some time t*. Such 
points, that we denote by Xq , are local minima of the ini- 
tial velocity gradient which have to be negative, so that 
the following relations are satisfied: 



d 2 * * 1 d 3 * 



dx 



o 



dx 3 



0, 



^(4)<0. (4.4) 



There is of course an extra global regularity condition 
that the preshock Lagrangian location Xq has not been 
captured by a mature shock at a time previous to t*. 
This global condition affects only constants but not the 
scaling behavior of p(p) at large p. 

We now Taylor-expand the initial density and the initial 
velocity potential in the vicinity of Xq. By adding a suit- 
able constant to the initial potential, shifting Xq to the 
origin and making a Galilean transformation canceling 
the initial velocity at Xq, we obtain the following "nor- 
mal forms" for the Lagrangian potential (3.5) and for 
the density 



®(xo,t)^Txl+Cx%, p(X(x ,t),t): 



-Po 



t+12(x^ 



(4.5) 



t-U 



and £ 



t* d^o 
24 dx* 



< 0. 



(4.6) 



Xq = 



The Lagrangian potential bifurcates from a situation 
where it has a single maximum at r < through a de- 
generate maximum with quartic behavior at r = 0, to 
a situation where convexity is lost and where it has two 
maxima at x = ±y / —r/ (4£) for r > 0. As a result of 
our choice of coordinates, the symmetry implies that the 
convex hull contains a horizontal segment joining these 
two maxima (see. figure 14(a)). 
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Fig. 14. Normal form of the Lagrangian potential, (a) in 
one dimension, in the time-neighborhood of a preshock; at 
the time of the preshock (t = 0), the Lagrangian potential 
changes from a single extremum to three extrema and devel- 
ops a non-trivial convex hull (shown as a dashed line), (b) 
in two dimension, the space neighborhood of a shock ending 
point has a structure similar to the spatio-temporal normal 
form of a preshock in one dimension when replacing the £0,2 
variable by the time r; the continuous line is the separatrix 
between the regular part and the ruled surface of the convex 
hull; the dotted line corresponds to the locations where the 
Jacobian of the Lagrangian map vanishes. 

We see from (4.5) that the Eulerian density p is propor- 
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tional to Xq in Lagrangian coordinates at t = t*. Since 
X = —d Xo ^, the relation between Lagrangian and Eu- 
lerian coordinates is cubic, so that at r = 0, the den- 
sity has a singularity oc x~ 2 / 3 in Eulerian coordinates. 
At any time t ^ t*, the density remains bounded ex- 
cept at the shock position. Before the preshock (r < 0), 
it is clear that p < —po/r, while after (r > 0), exclu- 
sion of the Lagrangian shock interval [xq ,Xq] implies 
that p < po/(2t). Clearly, large densities arc obtained 
only in the immediate space-time neighborhood of the 
preshock. More precisely, it follows from (4.5) that hav- 
ing p(x, t) > /i requires simultaneously 

|r| < ^ and \x\ < (-12C)- 1/2 (^) ^ (4.7) 

The tail of the cumulative probability of the density can 
be determined from the fraction of Eulerian space-time 
where p exceeds a given value. This leads to 

P>{p; x, t) = Prob{ ( o(x, t) > p] ~ C(x, t)p-^ 2 7 (4.8) 

where the constant C can be expressed as 

C(x,t)=At f |Cr 1/2 p 3 (M,C)dC, (4.9) 

J — oo 

A is a positive numerical constant and designates the 
joint probability distribution of the preshock space-time 
position and of its "strength" coefficient ( (see [54] for 
details). This algebraic law for the cumulative probabil- 
ity implies that the PDF of the mass density has a power- 
law tail with exponent —7/2 at large values. Actually 
this law was first proposed in [37] for the large-negative 
tail of velocity gradients in one-dimensional forced Burg- 
ers turbulence, a subject to which we shall come back in 
section 7. 

In higher dimensions it was shown in [54] that the main 
contribution to the probability distribution tail of the 
mass density does not stem from preshocks but from 
"kurtoparabolic" points. Such singularities (called ^3 
according to the classification of [62], which is summa- 
rized in section 2.3) correspond to locations which be- 
long to the regular part of the convex hull of the La- 
grangian potential $(x ,t) and where its Hessian van- 
ishes. The name kurtoparabolic comes from the Greek 
"kurtos" meaning "convex". These points are located 
on the spatial boundaries of shocks and generically form 
space-time manifolds of co-dimension 2 (persisting iso- 
lated points for d = 2, lines for d = 3, etc.). As in 
one dimension, the normal form of such singularities is 
obtained by Taylor-expanding in a suitable coordinate 
frame the Lagrangian potential to the relevant order 

$(f ,0-C4i+ E (-y^+^o,!^), (4.10) 

2<j<d 



where the different coefficients satisfy inequalities that 
ensure that the surface is below its tangent plane at 
Xq = 0. The typical shape of the Lagrangian potential in 
two dimensions is shown in figure 14(b). The positions 
where the Jacobian of the Lagrangian map vanishes can 
be easily determined from this normal form. The convex 
hull of $ and the area where the mass density exceeds the 
value fj, can also be constructed explicitly. An important 
observation is that, in any dimension, the scalar product 
of the vector {Jq = (xo.2, • • • , £o,d) with the vector (3 — 
(/?2, • • • , (3d) plays locally the same role as time does in 
the analysis of one-dimensional preshocks. 

When p — > 00, the cumulative probability can be esti- 
mated as 

P > (p; x,t) oc p~ 3/2 x p- 1 

from x ,i from 0-y o 

X 1 X ■ ■ ■ X 1 X ^1^. (4.11) 

from other components of y from time 

The only non-trivial contributions come from rro.i and 
from the component of yb along the direction of (3, all 
the other components and time contributing order-unity 
factors. Hence, the cumulative probability P > {p) is pro- 
portional to /i~ 5 / 2 in any dimension, so that the PDF of 
mass density has a power-law behavior with the univer- 
sal exponent —7/2. 

As we have seen, the theory is not very different in one 
and higher dimension even if kurtoparabolic points are 
persistent only in the latter case. This is due to the pres- 
ence of a time-like direction in the case d > 2. 

4-2 Evolution of matter inside shocks 

As we have seen in the previous subsection, the mass 
density becomes very large in the neighborhood of kur- 
toparabolic points (A 3 singularities) corresponding to 
the space-time boundaries of shocks. Such singularities 
dominate the tail of the mass density probability dis- 
tribution and contribute a power-law behavior with ex- 
ponent —7/2. However the mass distribution depends 
strongly on what happens inside the shocks where the 
density is infinite. Indeed, after the formation of the first 
singularity a finite fraction of the initial mass gets con- 
centrated inside these low-dimensional structures. De- 
scribing the mass distribution requires understanding 
how matter evolves once concentrated in the shocks. But 
before it will be useful to explain briefly the time evolu- 
tion of the shock manifold. 

4-2.1 Dynamics of singularities 

Suppose that X(t) denotes the position of a shock at 
time t. We suppose this singularity to be of type A\ 
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(see section 2.3), so that at this position, the veloc- 
ity field is discontinuous; we denote by v\ , . . . , v n the n 
different limiting values it takes at that point. At any 
time we generically have n < d + 1 and occasionally 
n = d + 2 at the space-time positions of shock metamor- 
phoses corresponding to instants when two Af singular- 
ities merge. We first restrict ourselves to persistent sin- 
gularities, meaning that n < d+ 1. In the neighborhood 
of X(t), it is easily checked that the velocity potential 
can be written as 



^(x,t) = $(X(t),t) + max Vj ■ (X(t) - x) 

j=~l. .n L 

+o(\\x-X(t)\\). 



(4.12) 



This expansion divides locally the physical space in n 
subdomains Clj where vj ■ (X(t) — x) is maximum, i.e. 



approximation of the potential at time t + St: 



+ max max 



v j .(X(t)-y)-—\\x-y\ 



(4.16) 



Note that here, St and \\x — X(t)\\ are chosen sufficiently 
small in a suitable way to ensure that (i) any singularity 
of higher co-dimension does not interfere with the posi- 
tion of X(t) between times t and t + St and that (ii) the 
subleading terms are always dominated by the kinetic 
energy contribution \\x — y|| 2 /(2<Ji). 

The two maxima in y and in j of (4.16) can be inter- 
changed, under the condition that the maximum in y is 
restricted to the domain flj defined in (4.13). The po- 
tential at time t + St can thus be written as 



(y-X(t)) > 0, 1 < i< n. (4.13) 



Writing the expansion (4.12) amounts to approximating 
the velocity potential by a continuous function which is 
piecewise linear on the subdomains flj . The boundaries 
between the fij's define the local shock manifold. The 
maximum in (4.12) ensures that we are focusing on en- 
tropic solutions to the Burgers equation (solutions ob- 
tained in the limit of vanishing viscosity) and results in 
the convexity of the local approximation of the poten- 
tial. Note also that the position x = X(t) of the reference 
singular point corresponds by construction to the unique 
intersection of all subdomains Qj. Remember that we 
have assumed that locally, the solution does not have 
higher-order singularity. 



+ ~ V(X(t),t) 



+ max max 



(X(t)-y)-l- t \\x-yr 



■ (4-17) 



We next remark that for all x, j and y, one has 



v r (X{t)-y) - — \\x-yf 



< Vj ■ (X(t) -x + StVj) 



St, 



u,|| 2 (,4.18) 



which gives an upper-bound to the maximum over y E 
Qj in (4.17). Suppose now that the maximum over the 
index j is achieved for j = j . This means that for all 
1 < i < n and y G fli 



The approximation (4.12) fully describes the local struc- 
ture of the singularity. If n = 2, corresponding to X{t) 
being the position of a simple shock, it is easily checked 
from (4.12) that there will actually exist a whole shock 
hyper-plane given by the set of positions y satisfying 



(v 1 -v 2 )-(X(t)-y)=0. 



(4.14) 



If n > 2, meaning that X(t) is an intersection be- 
tween different shocks, all the singular manifolds of 
co-dimension m < n are present in the expansion and 
are given by the set of positions y satisfying 



< max 

z6fij n 



rj j0 .(X(t)-z)-—\\x-z\\ 
St, 



<v jo -(X(t)-x + Stv jo )--\\v 



(4.19) 



Let Qi be the domain containing the vector (x — Stvj ). 
Then, (4.19) applied to i = io and y = x — Stvj trivially 
implies that 



(«io - «jo) ' (x-Stv jo -X(t))>0, 



(4.20) 



v n ■ (X(t) -y) = ...=v im . (X(t) - y) , (4.15) 
with 1 < i\ < ■ ■ ■ < i m < n. 



which together with the definition (4.13) for Q,i a leads 
to i{) = jo- Hence, to summarize, if the first maximum 
is reached for j = j then the second maximum is nec- 
essarily reached for y = x — Stvj . 



We next apply the variational principle (3.2) in order to 
solve the decaying problem between times t and t + St 
with the initial condition given by (4.12). This yields an 



It is clear that the approximation (4.16) of the velocity 
potential at time t + St preserves the local structure of 
the singular manifold. Indeed, for m < n, the positions 
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y satisfying 



St, 



v 1 .(X(t)-y)+"-\\v 1 f 



,v m -(X(t)-y) + j\\v m \\ 2 (4.21) 



form a (d — TO)-dimensional shock manifold. The trajec- 
tory X(t) of the reference singular point satisfies 



dX 1, 



Vl 



dt 2' 
which can be rewritten as 



dX 

~dT 



l, 



WdX/dt-v^l = ■■■ = \\dX/dt-V v 



(4.22) 



(4.23) 



This gives n relations for the d components of the vector 
dX/dt. These relations allow determining the normal 
velocity of the singular manifold. The tangent velocity 
remains undetermined. The velocity of the singularity 
located at X(t) is completely determined only if n = d, 
i.e. for point singularities. For instance when d = 1, the 
dynamics of shocks is given by 



dX 1 

^ = 2 (U1+U2) - 



(4.24) 



meaning that they move with a velocity equal to the half 
sum of their right and left velocities. For d — 2, only the 
positions of triple points (singularities of type A\ corre- 
sponding to the intersection of three shock lines) are well 
determined. It is easily checked that the two-dimensional 
velocity vector dX /dt is the circumcenter of the trian- 
gle formed by the three limiting values {v\, V2, V3) that 
are achieved by the velocity field at this position (see 
figure 15). 




Fig. 15. Determination of the velocity of a triple point and 
of that of the mass inside it when the three limiting values 
of the velocity Hi, V2, and S3 form an obtuse triangle. The 
dash-dotted circle is the circumcircle whose center gives the 
velocity of the singularity and the dashed circle is the small- 
est circle containing the triangle whose center gives the ve- 
locity of mass. 



4-2.2 Dynamics of the mass inside the singular mani- 
fold 

One of the central themes of this review article is a con- 
nection between Lagrangian particle dynamics and the 
inviscid Burgers equation. In the unforced case the ve- 
locity is conserved along particle trajectories minimizing 
the Lagrangian action (see section 2). At a given mo- 
ment of time, all particles whose trajectories are not min- 
imizers have been absorbed by the shocks. In the one- 
dimensional case when shocks are isolated points, par- 
ticles absorbed by shocks just follow the dynamics of a 
shock point. However, in the multi-dimensional case the 
geometry of the singular shock manifold can be rather 
complicated. This results in a non-trivial particle dy- 
namics inside the singular manifold. In other words, the 
particle absorbed by shocks have a rich afterlife and 
the main problem is to describe their dynamical proper- 
ties inside the singular manifold. This problem was ad- 
dressed by I. Bogaevsky in [18]. 

The basic idea is to consider first particle transport by 
the velocity field given by smooth solutions to the viscous 
Burgers equation. Indeed, let v"(x,t) be a solution to 
the viscous Burgers equation 

d t v" + (tT • V)tT = vV 2 v v - VF(x, t). 

Then the dynamics of a Lagrangian particle labeled by 
its position xq at time t — is described by the system 
of ordinary differential equations 

X v (x ,t) = v»{X v (x ,t),t), X v (x ,$)=x , (4.25) 

where the dots stand for time derivatives. It is possible to 
show that in the inviscid limit v — > solutions to (4.25) 
converge to limiting trajectories {X(xo,t}}. These lim- 
iting trajectories are not disjoint anymore. In fact, two 
trajectories corresponding to different initial positions 
Xq and Xq can merge: X(x\,t*) = X(xq, t*). This corre- 
sponds to absorption of particles by the shock manifold. 
Of course, two trajectories coincide after they merge: 
X(xq, t) = X(xq, t) for t >t*. Particles which until time 
t never merged with any other particles correspond to 
minimizers. Such trajectories obviously satisfy the lim- 
iting differential equation: 

X(x ,t)=v(X(x Q ,t),t), X(x Q ,0)=x , (4.26) 

where v(x, t) is the entropic solution of the inviscid Burg- 
ers equation which is well defined outside of the shock 
manifold. However, we are mostly interested in the dy- 
namics of particles which have merged with other par- 
ticles and thus were absorbed by shocks. One can prove 
that for such trajectories one-sided time derivatives exist 



km *<t 
dt At^0+ 



At) - X(t) 



At 



(4.27) 
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and satisfy a "one-sided" differential equation: 



Hence, the difference of actions can be written as 



(4.28) 



Here v^ s \-,t) is the velocity field on the shock manifold 
(index s stands for shocks) . ft turns out that v^ s ' 1 (x, t) 
and the corresponding shock trajectories satisfy a vari- 
ational principle, described hereafter. Denote by ^(x, t) 
a potential of the viscous velocity field v(x, t): v(x, t) = 
— V l &(a f , t). As we have pointed out many times be- 
fore, — ^(x, t) corresponds to a minimum Lagrangian ac- 
tion among all the Lagrangian trajectories which pass 
through point x at time t. Shocks correspond to a situa- 
tion where the minimum is attained for several different 
trajectories. Correspondingly, one has several smooth 
branches such that *&(x, t) = &i(x, t), 1 < i < k. Sup- 
pose a particle moves from a point of shock (x, t) with 
a velocity v. Then at infinitcsimally close time t + St its 
position will be x + vSt. In linear approximation (see pre- 
vious subsection) the Lagrangian action of this infinites- 
imal piece of trajectory is equal to [|w| 2 /2 — F(x,t)]St. 
Of course, the action minimizing trajectory at the point 
(x + vSt, t + St) does not pass through a shock point 
(x, t). Hence, the minimum action — ^(x + v5t,t + St) 



is smaller than — ^(x, t) 



2 /2 - F(x,t)]8t for any 



velocity v. However, we can put a variational condition 
on v which requires the difference between — *ff(x, t) + 
[||w|| 2 /2 - F(x,t)]St and -*(x + vSt,t + St) to be as 
small as possible. This is exactly the variational principle 
which determines the velocity v — (x, t) at a shock 
point. It is easy to see that in linear approximation 



AA = -V(x, t) + -\\v\\ 2 St + ^(x+vSt,t+St) 



= -\\v\\ 2 5t- min 

A Ki<k 



Vi-v- -\\vi\\ 2 



St 



= — max \\v — vA\ 2 St. 
2 i<i<fc 



(4.32) 



Obviously minimization of AA over v corresponds to a 
center of a minimum ball covering V{. It implies that 
such a center gives the velocity v^ s '(x(t), t) of particles 
concentrated at a shock point (x, t) . It is interesting that 
this variational principle implies that a particle absorbed 
by a shock cannot leave the singular shock manifold in 
the future. 

Let us now consider the first nontrivial generic exam- 
ple of a shock point, namely a triple point in two di- 
mensions d = 2. The point (X(t), t) is thus the intersec- 
tion of three shock lines. In this case there are exactly 
three smooth branches ^>i(-,t) with limiting velocities 
Vi = — V^j, 1 < i < 3. As we have seen in previous sec- 
tion the motion of the triple point is determined by con- 
tinuity of the velocity potential at (X, t). The "geomet- 
rical velocity" dX/dt of the triple point is then the cir- 
cumcenter of the triangle formed by the three velocities 
vi, V2, v 3 . It is easy to see that dX/dt — only in the 
case when the vectors v\, V2, and V3 form an acute trian- 
gle. If so, a cluster of particles follows the triple point. In 
the opposite case when the triangle is obtuse, the parti- 
cles leave the node. Such a situation is presented in fig- 
ure 15, where the mass leaves the node along the shock 
line delimiting the values v\ and V3 of the velocity. 



^(x+vSt,t+5t) = max{^i(x+v5t,t+6t)] 

= V(x,t)- min [-Wi{x,t) ■ v - d t Vi(x,t)}dt. (4.29) 

Ki<k 



Let us denote by Vi the limiting velocities — V*i(f, t) 
at the shock point (x, t) . Then, using Hamilton-Jacobi 
equation for the velocity potential 



d t <S> l (x,t) = ^\\W t (x,t)\\ 2 + F(x,t) 



— \\Vi 

2 II » 



+ F{x,t) 



(4.30) 



we have 



V(x + v6t, t + St) = t) 



mm 

Ki<k 



Vi-v- -\\vi\\ 2 



St-F(x, t)5t. (4.31) 



The analysis presented above was carried out for the 
Burgers equation jointly with A. Sobolevskh as a part of 
ongoing work on a similar theory for the case of a general 
Hamilton-Jacobi equation, with a Hamiltonian that is 
convex in the momentum variable. The formal extension 
of this analysis to the latter case is straightforward and 
can be left to the interested reader; however at present a 
rigorous justification of it, employing methods of [18], is 
known only for the case of H(x, x, t) — a(x, t)\x\ 2 , with 
a(x,t) > 0. 

4-3 Connections with convex optimization problems 

As discussed in section 1.2, Burgers dynamics is known 
in cosmology as the adhesion model and frequently used 
to understand the formation of the large-scale structures 
in the Universe. Recently, this model was used as a basis 
for developing new techniques for one of the most chal- 
lenging questions in modern cosmology, namely recon- 
struction. This problem aims at reconstructing the dy- 
namical history of the Universe through the mass den- 
sity initial fluctuations that evolved into the distribution 
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of matter and galaxies which is nowadays observed (see, 
e.g., [98]). The main difficulty encountered is that the 
velocities of galaxies (the peculiar velocities) are usually 
unknown, so that most approaches lead to non-unique 
solutions to this ill-posed problem. The reconstruction 
technique we present here, which was proposed in [55,25] , 
is based on the observation that, to the leading order, 
the mass is initially uniformly distributed in space (see, 
e.g., [98]). This observation, together with the Zeldovich 
approximation, leads to a reformulation of the prob- 
lem as a well-posed instance of an optimal mass trans- 
portation problem between the initial (uniform) and the 
present (observed) distributions of mass. More precisely 
it amounts to a convex optimization problem related 
to the Monge-Ampere equation and dually, as found 
by Kantorovich [73], to a linear programming problem. 
This is the reason why the name MAK (Monge-Ampere- 
Kantorovich) has been proposed for this method in [55]. 
Namely, one has to find the transformation from initial 
to current positions (the Lagrangian map) that maps 
the initial density p(x , 0) = po to the field p(x, t) which 
is nowadays observed. One then use a well-known fact 
in cosmology: because of the expansion of the Universe, 
the initial velocity field of the self-gravitating matter is 
slaved to the initial gravitational field (see, e.g., [25]). 
This observation implies that the initial velocity field is 
potential and allows one to deduce from it the sublead- 
ing fluctuations of the mass density. 

The MAK reconstruction technique is based on two cru- 
cial assumptions. First the Lagrangian map xo i— » x = 
X(xo, t) is assumed to be potential, i.e. X — V Xo <&(:ro). 
Second, the Lagrangian potential <&{xo) is assumed to 
be a convex function. As explained in [25] these two hy- 
potheses are motivated by the adhesion model (and thus 
inviscid Burgers dynamics) where they are trivially sat- 
isfied. As we will see later the reverse is actually true: the 
potentiality of the Lagrangian map and the convexity of 
the potential is equivalent to assuming that the latent ve- 
locity field is a solution to the Burgers equation. We will 
now see how, under these hypotheses, the reconstruction 
problem relates to Monge-Ampere equation. Conserva- 
tion of mass trivially implies that p(x, t)d 3 x — p d 3 xo, 
which can be rewritten in terms of the Jacobian matrix 
{dX^Kdxi) as 



dct 



dX l 
dx 3 



Po 



p(X{x ,t),t) 
Potentiality of the Lagrangian map leads to 
9 2 $ 



(4.33) 



det 



dx Q dx 3 



Po 



p(v* *,t)" 



(4.34) 



The problem with this formulation is that the unknown 
potential <& enters the right-hand side of the equation 



in a non-trivial way. Convexity of the Lagrangian po- 
tential <& is next used to reformulate the problem in 
term of the inverse Lagrangian map. Indeed, if $ is con- 
vex, the inverse Lagrangian map is also potential, i.e. 
xo = Xo(x,t) = \7 x Q(x) with the potential 9 itself con- 
vex. The two potentials $ and 9 are moreover related 
by Legendre transforms: 



Q(x) — max[f • So — $(a?o)], 

Xq 

$(xq) = max[x • xq — 9(x)]. 



(4.35) 
(4.36) 



In terms of the inverse Lagrangian potential 9 the con- 
servation of mass (4.34) reads 



dct 



d 2 e 



dx l dxi 



p(x,t), 



(4.37) 



which is exactly the elliptic Monge-Ampere equation. 
This time, the difficulty expressed above has disappeared 
since the unknown potential 9 does not enter the right- 
hand side of the equation. Note that we have implic- 
itly assumed here that the present distribution of mass 
has no singularity. The case of a singular distribution 
could actually be treated using a weak formulation of 
the Monge-Ampere equation, which amounts to apply- 
ing conservation of mass on any subdomain but requires 
allowing the inverse Lagrangian map to be multival- 
ued. The next step in the design of the MAK method is 
to reformulate (4.37) as an optimal transport problem 
with quadratic cost. Indeed, as shown in [24], the map 

X(xq, t) (and its inverse Ao(x, t)) minimizing the cost 



(4.38) 



T = J \\X(x ,t) - x \\ 2 p d 3 x 
= J \\x-X a (x, t)\\ 2 p(x, t) d 3 x, 



is a potential map whose potential is convex and is the 
solution to the Monge-Ampere equation (4.37). This can 
be understood using a variational approach as proposed 
in [55] . Suppose we perform a small displacement SXq (x) 
of the inverse Lagrangian map X {x,t) solution of the 
optimal transport problem. On the one hand the only ad- 
missible displacement are those satisfying the constraint 
to map the initial density field po to the final one p{x, t). 
It is shown in [25] that this is equivalent to require that 
\7 X ■ [p(x,t)5X (x)] = 0. On the other hand one easily 
see that the variation of the cost function corresponding 
to the variation Sx reads 



51= -2 



X (x, t)} ■ [p{x, t)SX (x)] d 3 x. (4.39) 



This integral can be interpreted as the scalar product 
(in the L 2 sense) between x — A (x , t) and p(x)SX(xo)- 
Hence the optimal solution, which should satisfy SI = 
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for all SXq, is such that the displacement x — Xo(xo,t) 
(or equivalently X(x ) — x*o) is orthogonal to all 
divergence-free vector fields. This means that it is nec- 
essarily the gradient of a potential, from which it follows 
that X(xo,t) = V Xo Q(xo). Convexity follows from the 

observation that the Lagrangian map Xo i— > X has to 
satisfy 



(f - 4) ' [*(*<>) - Xfo)] > 0. 



(4.40) 



Indeed, if that was not the case, one can easily check that 
any map where the Lagrangian pre- image of a neighbor- 
hood of x*o and of one of x*o are inverted would lead to a 
smaller cost. Formulated in terms of potential maps, the 
relation (4.40) straightforwardly implies convexity of $. 
This finishes the proof of equivalence between Monge- 
Ampere equation and the optimal transport problem 
with quadratic cost. 

The goal of reformulating reconstruction as an optimiza- 
tion problem is mostly algorithmic. Once discretized, 
the problem of finding the optimal map between initial 
and final positions amounts is equivalent to solving a so- 
called assignment problem. An efficient method to deal 
numerically with such problems is based on the auction 
algorithm [15] and was used in [25] with data stemming 
from iV-body cosmological simulations. As summarized 
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Fig. 16. Test of the MAK reconstruction for a sample of 
N' = 17, 178 points from a iV-body simulation (from [25]). 
The scatter diagram plots reconstructed versus true initial 
positions. The histogram inset gives the distribution (in per- 
centages) of distances between true and reconstructed initial 
positions; the horizontal unit is the distance between two 
sampled points. The width of the first bin is less than unity to 
ensure that only exactly reconstructed points fall in it. More 
than sixty percent of the points are exactly reconstructed. 

in figure 16, the MAK reconstruction method leads to 
very promising results. More than 60% of the discrete 
points are assigned to their actual Lagrangian pre-image. 
Such a number has to be compared with other recon- 



struction methods for which the success rate barely ex- 
ceed 40% for the same data set. 

Even if the mapping from initial to final positions is 
unique, the peculiar velocities are not well defined ex- 
cept if we have some extra knowledge of what is happen- 
ing at intermediate times < t' < t. Of course the den- 
sity field p(x',t') is unknown. However, there are triv- 
ial physical requirements. First the two mass transport 
problems between and if and between if and t have 
both to be optimal. This means that one looks for two 
Lagrangian maps, X\ from to t' and X 2 from t' to t 
which are minimizing the respective costs 



= J \\Xi(x ) -x \\ 2 p d 3 x , 
= J \\x-X^{x)\\ 2 p{x,t) d 3 x. 



(4.41) 



The second physical requirement is that the composition 
of these two optimal maps have to give the Lagrangian 
map between times and t, namely X(x , t) = A 2 (Ai). 
Under these two conditions there is equivalence between 
the optimal transport with a quadratic cost and the 
Burgers dynamics supplemented by the transport of a 
density field (see [13] for details). 

5 Forced Burgers turbulence 

5.1 Stationary regime and global minimizer 

We consider in this section solutions to the forced Burg- 
ers equation. As we have seen in section 2, the solution in 
the limit of vanishing viscosity can be expressed at any 
time t in terms of the initial condition at time to through 
a variational principle which consists in minimizing an 
action along particle trajectories. The statistically sta- 
tionary regime toward which the solution converges at 
large time can be studied assuming that the by reject- 
ing the initial time to is at minus infinity. The solution 
is then given by the variational principle 



t) = 



-inf 

7(0 



da}, (5.1) 



where the infimum is taken over all (absolutely continu- 
ous) curves 7 : (—00, t] — ► fl such that j(t) = x. In this 
setting, the action is computed over the whole half line 
(— 00, t] and the argument of the infimum does not de- 
pend anymore on the initial condition. Of course, (5.1) 
defines VP up to an additive constant. This means that 
only the differences ^(x, t) — ^(0, t) can actually be de- 
fined. A trajectory 7 minimizing (5.1) is called a one- 
sided minimizer. It is easily seen from (5.1) that all the 
minimizers are solutions of the Euler-Lagrange equation 



7>) = -VF(7(«),«) 



(5.2) 
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where the dots denote time derivatives. This equation 
defines a 2<i-dimensional (possibly random) dynamical 

system in the position- velocity phase space (7, 7) . The 
Lagrangian one-sided minimizers 7 defined over the half- 
infinite interval (— 00, t] play a crucial role in the con- 
struction of the global solution and of the stationary 
regime. Namely, a global solution to the randomly forced 
inviscid Burgers equation is given by v(x, t) = j(t) where 
j(t) = x. To prove that such half-infinite minimizers ex- 
ist, one has to take the limit to — > —00 for minimizers 
defined on the finite time interval [to,t]. The existence of 
this limit follows from a uniform bound on the absolute 
value of the velocity |t*| (see, e.g., [38]). Obtaining such 
a bound becomes the central problem for the theory, as 
we shall now see. 

When the configuration space f2 where the solutions live 
is compact (bounded), one can expect the velocity of a 
minimizer to be uniformly bounded. Indeed, in this case 
the displacement of a minimizer for any time interval is 
then bounded by the diameter of the domain Q, so that 
action minimizing trajectories cannot have large veloc- 
ities. For forcing potential that are delta-correlated in 
time, it has been shown by E et al. [38] in one dimension 
and by Iturriaga and Khanin [68,69] in higher dimen- 
sions that the minimizing problem (5.1) has a unique 
solution \I> with the following properties: 

• ^ is the unique statistically stationary solution to the 
Hamilton-Jacobi equation (2.2) in the inviscid limit 
v -> 0; 

• is almost everywhere differentiable with respect to 
the space variable x; 

• — uniquely defines a statistically stationary solu- 
tion to the Burgers equation in the inviscid limit; 

• there exists a unique one-sided minimizer at those Eu- 
lerian positions x where the potential \& is differen- 
tiable; the locations where ^ is not differentiable cor- 
respond to shocks. 

• There exists a unique minimizer 7^ that minimizes 
the action calculated from —00 to any time t. It is 
called the global minimizer (or two-sided minimizer) 
and corresponds to the trajectory of a fluid particle 
that is never absorbed by shocks. Moreover, all one- 
sided minimizers are asymptotic to it as s — > —00. 

All the properties above follow from the variational ap- 
proach. In fact, the variational principle (2.12) imply 
similar statements in the viscous case. Of course, when 
viscosity is positive the unique statistically stationary 
solution is smooth. However, one can show that the sta- 
tionary distribution corresponding to such solutions con- 
verges to inviscid stationary distribution in the limit 
v — > [58] . Although the variational proofs are concep- 
tual, general and simple, they are based on the fluctua- 
tion mechanism and therefore do not give a good control 
of the rate of convergence to the statistically stationary 
regime. Exponential convergence would follow from the 



hyperbolicity of the global minimizer. Although one ex- 
pects hyperbolicity holds in any dimension, mathemat- 
ically it is an open problem. At present a rigorous proof 
of hyperbolicity is only available in dimension one [38] . 

The assumption of compactness of the configuration 
space f2 is essential in the construction of the stationary 
regime. As we will see in subsection 5.4, the situation 
is much more complex in the non-compact case when 
for instance the solution is defined on the whole space 
fi = R d . 

5.2 Topological shocks 

To introduce the notion of topological shock we first fo- 
cus on the one-dimensional case in a periodic domain, 
i.e. in O = T = R/Z. If we "unwrap" at a given time t 
the configuration space to its universal cover R (see fig- 
ure 17(a)), we then obtain an infinite number of global 

minimizer 7^ , which at all time s < t satisfy 7^ (s) = 

jj^\s) + 1. All the one-sided minimizers converge back- 
ward in time to one of these global minimizers. The topo- 
logical shock (or main shock) is defined as the set of x 
positions giving rise to several minimizers approaching 
two successive replicas of the global minimizer. This par- 
ticular shock is also the only shock that has existed for 
all times. 

This construction can easily be extended to higher 
dimensions (sec [10]). For this we unwrap the d- 
dimensional torus T d to its universal cover, the full 
space R d (see figure 17(b) for d = 2). Then, the different 
replicas of the periodic domain define a lattice of global 
minimizers parameterized by integer vectors k. The 
backward-in-time convergence on the torus of the one- 
sided minimizers to the global minimizer implies that 
a minimizer associated to a location x in R d at time t 
will be asymptotic to one of the global minimizer of 
the lattice. Hence, every position x which has a unique 
one-sided minimizer is associated to an integer vector 
k(x). This defines a tiling of space at time t. The tiles 
Oj: are the sets of points whose associated one-sided 

minimizers are asymptotic to the k-th global minimizer. 
The boundaries of the O^'s are the topological shocks. 
They are the locations from which at least two one- 
sided minimizers approach different global minimizers 
on the lattice. Indeed, a point where two tiles 0% and 
0£ 2 meet, has at least two one-sided minimizers, one 

of which is asymptotic to ji^ and another to 7^' . Of 

fel K2 

course, there are also points on the boundaries where 
three or more tiles meet and thus where more than two 
one-sided minimizcrs are asymptotic to different global 
minimizers. For d — 2 such locations are gencrically iso- 
lated points corresponding to the intersections of three 
or more topological shock lines, while for d = 3, they 
form edges and vertices where shock surfaces meet. Note 
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Fig. 17. Space-time sketch of the unwraping of the periodic 
domain T to the whole space R d for d = 1 (a) and d — 2 (b). 

that, generically, there exist other points inside 0£ with 
several minimizers. They correspond to shocks of "lo- 
cal" nature because at these locations, all the one-sided 
minimizers are asymptotic to the same global minimizer 
7^ and hence, to each other. In terms of Lagrangian 
dynamics, the topological shocks play a role dual to that 
of the global minimizer. Indeed, all the fluid particles 
are converging backward-in-time to the global mini- 
mizer and are absorbed forward-in-time by the topo- 
logical shocks. For the transportation of mass when we 
assume that the Burgers equation is supplemented by 
a continuity equation for the mass density, all the mass 
concentrate at large times in the topological shocks. 

The global structure of the topological shocks is related 
to the various singularities generically present in the so- 
lution to the Burgers equation that were detailed in sec- 
tion 2.3. Generically there are no locations associated to 
more than (d+ 1) minimizers. As one expects to see only 
generic behavior in a random situation, the probability 
to have points with more than (d + 1) one-sided mini- 
mizers is zero. It follows that there are no points where 
(d + 2) tiles 0g meet, which is an important restriction 
on the structure of the tiling. For d = 2 it implies that 
the tiling is constituted of curvilinear hexagons. Indeed, 
suppose each tile 0£ is a curvilinear polygon with s ver- 
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Fig. 18. (a) Position of the topological shock on the torus; 
the two triple points are represented as dots, (b) Snapshot 
of the velocity potential ip(x,y,t) for d — 2 in the statisti- 
cal steady state, obtained numerically with 256 2 grid points. 
Shock lines, corresponding to locations where tp is not dif- 
ferentiable, are represented as black lines on the bottom of 
the picture; the four gray areas are different tiles separated 
by the topological shocks; the other lines are local shocks. 

tices corresponding to triple points. For a large piece 
of the tiling that consists of N tiles, the total num- 
ber of vertices is n v <~ sN/3 and the total number of 
edges is n e ~ sN/2. The Euler formula implies that 
1 = n v — n e + N ~ (6 — s)N/6, and we necessarily have 
s = 6, corresponding to an hexagonal tiling. As shown 
in figure 18(a), this structure corresponds on the peri- 
odicity torus T 2 , to two triple points connected by three 
shock lines that are the curvilinear edges of the hexagon 
Og. The connection between the steady-state potential 
and the topological shocks is illustrated numerically on 
figure 18(b). The different tiles covering the periodic do- 
main were obtained by tracking backward in time fluid 
particle trajectories and by determining to which peri- 
odic image of the global minimizer they converge. 

In dimensions higher than two, the structure of topolog- 
ical shocks is more complicated. For instance it is not 
possible to determine in a unique manner the shape of 
the polyhedra forming the tiling. However, it has been 
shown by Matveev [87] that for d = 3 the minimum poly- 
hedra forming such tiling has 24 vertices and 36 edges 
and is composed of 8 hexagons and 6 rectangles (see fig- 
ure 19). It is of interest to note that the structure of 
topological shocks is in direct relation with the notions 
of complexity and minimum spines of manifolds intro- 
duced by Matveev from a purely topological viewpoint. 
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Fig. 19. Sketch of the simplest configuration of the topolog- 
ical shock in dimension d = 3. 

Algebraic characterization of the topological shock 

In two dimensions, when periodic boundary conditions 
are considered, very strong constraints are imposed on 
the structure of the solution. In particular, the topol- 
ogy of the torus T 2 imply that the topological shocks 
generically form a periodic tiling of R 2 with curvilinear 
hexagons. However, this tiling can be of various alge- 
braic types. Consider the tile Ojj surrounded by its six 
immediate neighbors Oj: , where the integer vectors ki 

are labeled in anti-clockwise order, k\ having the small- 
est polar angle (see figure 20). It is easily seen that the 
periodicity of the tiling implies 

= ki~ ki, ki = —ki, h = -k 2 

and fc 6 = fci — fc 2 , (5.3) 

so that the whole information on the algebraic struc- 
ture of the tiling is contained in the vectors k\ and k 2 
which form a matrix S from the group SL(2, Z) of 2 x 2 
integer matrices with unit determinant. The matrix S 
gives information on the number of times each shock line 
turns around the torus before reconnecting to another 
triple point. Figure 18(a) corresponds to the simplest 
case when S is the identity matrix. When the forcing is 
stochastic, the matrix S is random and stationary solu- 
tions to the two-dimensional Burgers equation define a 
stationary distribution on SL(2, Z). 




Fig. 20. The algebraic structure of the topological shock in 
dimension d — 2 is determined by the indexes corresponding 
to immediate neighbors of the tiling considered. 



Certainly, topological shocks evolve in time and may 
change their algebraic structure. This happens through 
bifurcations (or metamorphoses) described in section 
2.3. In two dimensions, the generic mechanism which 
transforms the algebraic structure of topological shocks 
is the merger of two triple points. This metamorphosis 
is called the flipping bifurcation and corresponds to the 
appearance at time t* of an A\ singularity in the solu- 
tion associated to a position with four minimizers. The 
mechanism transforming the algebraic structure of the 
topological shock is illustrated in figure 21. Issues such 
as the minimum number of flips needed to transform the 
matrix S\ associated to the algebraic structure of the 
topological shock to another matrix 52 are discussed in 
in [1]. 



t<t t t = t t t>t. 




Fig. 21. Sketch of the tiling before, at the flipping time t* and 
after it. This example corresponds to a bifurcation from the 
matrix Si = [J i] to £2 = [1 "2]- The dashed boxes represent 
the periodicity domain [0, l] 2 . 

5.3 Hyperbolicity of the global minimizer 

The nature of the convergence to a statistical steady 
state is determined by the local properties of the global 
minimizer. The hyperbolicity of this action-minimizing 
trajectory implies an exponential convergence, so that 
the global picture of the solution is reached very rapidly, 
after just a few turnover times. 

Since the trajectory of the global minimizer is unique 
and can be extended to arbitrary large times, it corre- 
sponds to an ergodic invariant measure for the stochastic 
flow defined by the Euler-Lagrange equation (5.2). Con- 
ditioned by the random force, this measure is simply the 
delta measure sitting at the location (7^(0), 7^(0)). 
By the Osclcdcts ergodic theorem (see, e.g. [98]), 2c? non- 
random Lyapunov exponents can be associated to the 
global minimizer trajectory. Since the flow is symplectic 
these non-random exponents come in pairs with oppo- 
site signs. That is 

Ai>--->A d >0>-A d >--->-Ai. (5.4) 

Hyperbolicity is defined as the non- vanishing of all these 
exponents. Thus, the issue of hyperbolicity can be ad- 
dressed in terms of the backward-in-time convergence of 
the one-sided minimizers to the global one or, better, in 
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terms of forward-in-time dynamics. In the latter case, 
this amounts to looking how fast Lagrangian fluid par- 
ticles are absorbed by shocks. For this we consider the 
set fi reg (T) of locations x such that the fluid particle 
emanating from x at time t = survives, i.e. is not ab- 
sorbed by any shock, until the time t = T. The long-time 
shrinking of il rog as a function of time is asymptotically 
governed by the Lyapunov exponents. To ensure the ab- 
sence of vanishing Lyapunov exponents, it is sufficient to 
show that the diameter of f2 reg (T) decays exponentially 
as T — > +00. 

In one dimension, it has been shown in [38] that this is 
indeed the case, and particularly that there exists posi- 
tive constants a, /3, A and B such that 

Prob {diam f7 rcg (T) > A C - aT ) < Bc"^ . (5.5) 

Unfortunately this proof of hypcrbolicity is purely one- 
dimensional and at present time there is no extension of 
this result to higher dimensions. 

In two dimensions, the behavior of diam Q veg (T) at large 
times was studied numerically in [10] by using the fast 
Legendre transform described in section 2.4 and a forc- 
ing that is a sum of independent random impulses con- 
centrated at discrete times. The ideas of this numerical 
method are related to the Lagrangian structure of the 
flow. This easily permits to track numerically the set f2 rog 
of regular Lagrangian locations. As seen from figure 22, 
the diameter of this set decays exponentially fast in time 
for three different types of forcing, providing good ev- 
idence of the hyperbolicity of the global minimizer for 
d= 2. 
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Fig. 22. Time evolution of the diameter of the Lagrangian 
set fi(T) (points corresponding to the regular region) for 
three different types of forcing spectrum; average over 100 
realizations and with 256 2 grid points (from [10]). 



Hyperbolicity of the global minimizer implies existence 
at any time t of two d-dimensional smooth manifolds 
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Fig. 23. Sketch of the unstable manifold for d = 1 in the 
(x,v) plane. Shock locations (Af singularities) are obtained 
by applying Maxwell rules to the loops. A preshock (A3 
singularity) is represented; it corresponds to the formation 
of a loop in the manifold. The velocity profile which is the 
actual solution to the Burgers equation is represented as a 
bold line. 

in phase space (7,7) that are invariant by the Euler- 
Lagrange dynamics (5.2): a stable (attracting) manifold 
r( s \t) and an unstable (repelling) manifold r( u )(i), de- 
fined as the instantaneous location of trajectories con- 
verging to the global minimizer forward in time and 
backward in time, respectively. Since all the minimizers 
converge backward in time to the global minimizer, the 
graph in the position- velocity phase space (x, v) of the 
solution in the statistical steady state is made of pieces of 
the unstable manifold r( u ) (t) with discontinuities along 
the shocks lines or surfaces. In other words, shocks ap- 
pear as jumps between two different folds of the unstable 
manifold. The smoothness of the unstable manifold is an 
important property; for instance, it implies that when 
d = 2, the topological shock lines are smooth curves. 

In one dimension, where hyperbolicity is ensured, the 
main shock corresponds to a jump between the right 
branch and the left branch of the unstable manifold. Its 
position can be obtained geometrically after observing 
that the area b covered by the unstable manifold, once 
the latter is cut by the main shock, should be equal to 
the first integral of motion which is conserved, i.e. 

6 = y v(x,t)dx = J v {x)dx. (5.6) 

The other shocks (or secondary shocks that have existed 
only for a finite time) cut through the double-fold loops 
of the unstable manifold (see figure 23). Their locations 
can be obtained by a Maxwell rule applied to those loops. 
Indeed, the difference of the two areas defined by cutting 
such a loop at some position x is equal to the difference 
of actions of the two trajectories emanating from the up- 
per and lower locations and, thus, vanishes at the shock 
location. We will see in section 6 that this construction 
of the solution is also valid when the forcing is periodic 
in time, problem which can be related to Aubry-Mather 
theory relative to commensurate-incommensurate phase 
transitions. 
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The above geometrical construction of the solution has 
much in common with that appearing in the unforced 
problem. Indeed, as we have seen in section 3.1, when 
F = the solution can be obtained geometrically by 
considering in the (x, v) space, the Lagrangian manifold 
defined by the position and the velocity of the fluid par- 
ticles at a given time. This analogy gives good ground 
predicting that some universal properties associated to 
the unforced problem will still hold in the forced case?, as 
we will indeed see in section 7. Another instance concerns 
transport of mass in higher dimension. We have seen in 
section 4.1 that, for the unforced case, large but finite 
mass densities are localized near boundaries of shocks 
( "kurtoparabolic" singularities) contributing power-law 
tails with the exponent —7/2 to the probability density 
function of the mass density. When a force is applied the 
smoothness of the unstable manifold associated to the 
global minimizer should lead to the same universal law. 

5.4 The case of extended systems 

So far, we have discussed the global structure of the 
solution to the forced Burgers equation with periodic 
boundary conditions. Is is however of physical interest 
to understand instances when the size of the domain is 
much larger than the typical length scale of the forcing. 
In this section, we will focus on describing, in the one- 
dimensional case, the singular structure of the solution 
in unbounded domains. Based on the formalism of [11], 
we achieve this goal by considering a spatially periodic 
forcing with a characteristic scale L f much smaller than 
the system size L. More precisely, for a fixed size L we 
consider the stationary regime corresponding to the limit 
t — > 00 and then study the limit L — > 00 by keeping 
constant the energy injection rate (i.e. the C 2 norm of 
the forcing grows like L). 

In order to get an idea of the behavior of the solution, 
the limit of infinite aspect ration L/Lf was investigated 
numerically in [11]. As seen from figure 24(a) numerical 
observations suggest that at any time in the statistical 
steady state, the shape of the velocity profile is simi- 
lar to the order-unity aspect ratio problem, duplicated 
over independent intervals of size Lf. In particular, when 
tracking backward in time the trajectories of fluid par- 
ticles the minimizers converge to each other in a very 
non-uniform way. Figure 24(b) shows that the minimiz- 
ers form different branches, which are converging to each 
other backward in time; in space time a tree structure is 
obtained. As shown in figure 25(a) a similar behavior is 
observed for shocks. 

The velocity field at a given time t, consists of smooth 
pieces separated by shocks. Let us denote by {f2j} the 
set of intervals in [0, L), on which the solution u(-,t) 
is smooth. The boundaries of the Cl/s are the shocks 
positions. Each of these shocks is associated to a root- 
like structure formed by the trajectories of the various 
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Fig. 24. (a) Upper: snapshot of the velocity field for 
L = 256L/. Lower: zoom of the field in a interval of 
length 101//. (b) Minimizing trajectories in space time for 
L — 256 Lf and over a time interval of length T = 100 

shocks that have merged at times less than t to form the 
shock under consideration (see figure 25(a)). This root- 
like structure contains the whole history of the shock 
and in particular its age (i.e. the length of the deeper 
branch of the root structure). Indeed, if the root has a 
finite depth, the shock considered has only existed for a 
finite time. A T -global shock is defined as a shock whose 
associated root is deeper than — T. They can alterna- 
tively be defined geometrically by considering the left- 
most and the rightmost minimizer associated to it. After 
tracing them backward for a sufficiently long time, these 
two minimizers are getting close and eventually converge 
to each other exponentially fast (see figure 24(b)). For 
a T-global shock, the time when the two minimizers are 
getting within a distance smaller than the forcing cor- 
relation length Lf is larger than T. As we have seen in 
section 5.2, the existence in one dimension of a main 
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Fig. 25. (a) Shock trajectories for aspect ratio L/Lf = 32 
and with T = 10. The different gray areas correspond to 
the space-time domains associated to the different smooth 
pieces Qj of the velocity field at time t = 0. (b) Sketch of 
the space-time evolution of a given smooth piece Qj located 
between two shock trajectories X\ (t) and X2 (t) that merge 
at time Tj . 

shock in the spatially periodic situation follows from a 
simple topological argument. The main shock can also 
be defined as the only shock that has existed forever in 
the past. It is hence infinitely old, contrary to all other 
shocks, all of them being created at a finite time and 
having a finite age. When the periodicity condition is 
dropped, the main shock disappears and it is useful to 
consider the T-global shocks that mimic the behavior of 
a main shock over time scales larger than T. 

One can dually define T-global minimizers. All the 
smoothness intervals Slj defined above, except that 
which contains the global minimizcr, will be entirely 
absorbed by shocks after a sufficient time (sec figure 
25(b)). For each of these pieces, one can define a life-time 
Tj as the time when the last fluid particle contained in 
this piece at time t enters a shock. It corresponds to the 
first time for which the shock located on the left of this 
smooth interval at time t merges with the shock located 
on the right. When the life-time of such an interval is 



greater than T, the trajectory of the last surviving fluid 
particle is here called a T-global minimizcr. Note that, 
when T — > 00, the number of T-main shocks and of 
T-global minimizers is one, recovering respectively the 
notions of main-shock and of two-sided minimizer. 
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Fig. 26. Density of T-main shocks as a function of T for three 
different system sizes L/Lf = 64, 128 and 256; average over 
100 realizations. Lower inset: local scaling exponent. 

Hence, at a given instant t, and for any timelag T, 
the spatial domain [0, L) contains a certain number of 
T-objects. We define their spatial density as being the 
number of such objects, averaged with respect to the 
forcing realizations, divided by the size of the domain L. 
The density p(T) of T-main shocks was investigated nu- 
merically in [11] for the kicked case by using a two-step 
method: first, the simulation was run until a large time t 
for which the statistically stationary regime is reached; 
secondly, each shock present at time t was tracked 
backward-in-time down to the instant of its creation, 
giving an easy way to characterize the density p(T). It 
is seen in figure 26 that, for three different aspect ratios 
L/Lf, the density p(T) displays a power-law behavior 
p(T) oc T~ 2 / 3 for the intermediate time asymptotics 

rms- 

We now present a simple phenomenological theory aim- 
ing to explain the scaling exponent 2/3. We consider the 
solution at a fixed time (t = 0, for instance). Denote by 
£(T) the typical spatial separation scale for two nearest 
T-global shocks. Obviously, £{T) ~ l/p(T). The mean 
velocity of the spatial segment of length I is given by 

h = \ f u(x,0)dx (5.7) 

Since the expected value (u(x, 0)) = 0, and that the 
integral in (5.7) is over an interval of size much larger 
than the forcing correlation length, it is equivalent to 
a sum of independent centered random variables and 
scales as the Brownian motion. Hence, for large t one 
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has the following asymptotics 



than the forcing scale (see figure 27) . In order to observe 



/ u(x, 0) da; - VI, 
J[y,v+t] 



(5i 



which gives b e ~ £~ x l 2 for mean velocity fluctuations. 
Consider now the rightmost minimizer corresponding to 
the left T-global shock and the leftmost minimizer re- 
lated to the right one. Since there are no T-global shocks 
in between, it follows that the two minimizcrs we se- 
lected get close to each other backward-in-time around 
times of the order of — T. This means that the backward- 
in-time displacement of a spatial segment of length 0(£) 
is itself 0{£) for time intervals of the order of T. The 
corresponding displacement is given as the sum of two 
competing behaviors: the first, which can be understood 
as a drift induced by the local mean velocity bi, is due 
to the mean velocity fluctuations and is responsible for 
a displacement oc b{T; the second contribution is due to 
a standard diffusive scale oc T 1 / 2 expressing the diffu- 
sive behavior of the minimizing trajectories. Taking into 
account both terms we obtain 



£~ B 1 T£~ 1 / 2 + B 2 T 1 / 2 , 



(5.9) 



where B\ and B 2 are numerical constants. It is easy to 
see that the dominant contribution comes from the first 
term. Indeed, if the second term were to dominate, then 
£ would be much larger than T, which contradicts (5.9). 
Hence, one has £ <~ B{T£~ X I 2 , leading to the scaling 
behavior 



£{T) oc T 2/3 , p(T) oc T~ 



2/3 



(5.10) 



As we have already discussed, T-global shocks are shocks 
older than T. Denote by p(A) the probability density 
function (PDF) for the age of shocks. More precisely, 
p(A) is a density in the stationary regime of a probability 
distribution of the age A(t) of a shock, say the nearest 
to the origin. It follows from (5.10) that the probability 
of shocks whose age is larger than A decays like A~ 2 > 3 ; 
this implies the following asymptotics for the PDF p{A) : 



p{A) oc A" 5 / 3 . 



(5.11) 



Actually, the power-law behavior of the density p{T) of 
T-global shocks can be interpreted in term of an inverse 
cascade in the spectrum of the solution (although there 
is no conserved energy-like quantity). Indeed, the fluctu- 
ations (5.8) of the mean velocity suggest that, for large- 
enough separations I, the velocity potential increment 
scales like 



(x + e, t)-ip(x, t)\ ot£ 1/2 . 



(5.12) 



This behavior is responsible for the presence of an inter- 
mediate power-law range with exponent —2 in the spec- 
trum of the velocity potential at wavenumbers smaller 




Fig. 27. Spectrum {4> 2 (k)) of the velocity potential in the 
stationary regime for the aspect ratio L/Lf = 128. This 
spectrum contains two power-law ranges: at wavenumbers 
k ^> L/Lf, the traditional oc fc~ 4 inertial range connected 
to the presence of shocks in the solution and, for k <C L/Lf, 
an "inverse cascade" oc k~ 2 associated to the large-scale 
fluctuations of ip 

the k~ 2 range at small wavenumbers, the spectrum of 
the forcing potential must decay faster than k~ 2 ; other- 
wise the leading behavior is non-universal but depends 
on the functional form of the forcing correlation. 

The one-dimensional randomly forced Burgers equation 
in an unbounded domain has been studied in [66] with 
a different type of forcing: it was assumed that the forc- 
ing potential has at any time its global maximum and 
its global minimum in a prescribed compact region of 
space. It was proven that with these particular settings 
the statistically stationary regime exists and is very sim- 
ilar to that arising in compact domains. In particular, 
there exists a unique global minimizer located in a finite 
spatial interval for all times and all other minimizers are 
asymptotic to it in the limit t — > — oo. The main idea 
behind considering such type of forcing potential is to 
ensure that the potential energy plays a dominant role 
in comparison with the kinetic (elastic) term in the ac- 
tion. This leads to effective compactification and allows 
estimates on the velocities of fluid particles. As we al- 
ready mentioned in section 5.1, these estimates are very 
important and pave the way to the construction of the 
whole theory of the statistically stationary regime. 

Note finally that it was shown in [76] that for special 
cases of forcing potentials F(x, t), the velocity of a min- 
imizers can be arbitrarily large. More specifically, one 
can construct pathological forcing potentials such that 
minimizers are accelerated and reach infinite velocities. 
Randomness is of course expected to prevent such a type 
of non-generic blow-up. 
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6 Time-periodic forcing 



two times t > t' between which no force is applied: 



This section is devoted to the study of the solutions 
to the one-dimensional Burgers equation with time- 
periodic forcing. In this case many of the objects we 
have discussed above can be constructed almost explic- 
itly: the global minimizer, the main shock etc. Also, 
a mathematical analysis is then much simpler. For in- 
stance, hyperbolicity of the global minimizer follows 
immediately from first principles. Finally, the case of 
time-periodic forcing is directly related to the Aubry- 
Mather theory as we explain below. 

6. 1 Kicked Burgers turbulence 

We shall be concerned here with the initial-value prob- 
lem for the one-dimensional Burgers equation when the 
force is concentrated in Dirac delta functions at discrete 
times: 



f(x,t) 



]fj(x)6(t-tj), 



(6.1) 



ip(x,t) 



mm 

v 



(x - y) 2 

2{t-t>) 



4>(y,t') 



(6.6) 



As before, when t' is the initial time, the position y which 
minimizes (6.6) is the Lagrangian coordinate associated 
to the Eulerian coordinate x. The map y i— » x is called 
the Lagrangian map. By expanding the quadratic term it 
is easily shown that the calculation ofijj(-,t) from ip(-, t') 
is equivalent to a Legendre transformation. For details, 
see [104,107]. 

We now turn to the forced case with impulses applied at 
the kicking times tj. Let tj( t ) be the last such time before 
t. Using (6.6) iteratively between kicks and changing 
the potential ip(y,tj+i) discontinuously by the amount 
Fj + i{y) at times tj+i, we obtain 



ip(x,t) 



min [A{{yj};x,t;jo))-ip (y jo )], (6.7) 



{Vj}j <3< 



where both the "impulses" fj (x) and the "kicking times" 
tj are prescribed (deterministic or random) . The kicking 
times are ordered and form a finite or infinite sequence. 
The impulses fj(x) are always taken spatially smooth, 
i.e. acting only at large scales. The general scheme we 
are presenting below holds for any sequence of impulses 
fj(x) and kicking time. Later on we shall assume that 
they define a time-periodic forcing. The precise meaning 
we ascribe to the Burgers equation with such forcing is 
that at time tj, the solution u(x,t) changes discontinu- 
ously by the amount fj (x) 



u{x,t j+ ) = u(x,tj-) + fj(x), 



(6.2) 



while, between tj + and the solution evolves ac- 

cording to the unforced Burgers equation. 

We shall also make use of the formulation in terms of the 
velocity potential ip(x, t) and the force potentials Fj(x) 



u(x,t) = -d x iP(x,t), fj(x) = -—Fjix). (6.3) 
The velocity potential satisfies 

M = \{d x ^f + vd xx ^ + Fj(x) S(t -tj), (6.4) 



ip(x,t ) = ^o{x), 



(6.5) 



A({yj};x,t;j ) 



_ (x - yjf 



2(t - tj) 



.7-1 

+ E 

j=3o 



(y j+1 -yjf 



2(t j+1 - tj) 



- F 3+i(y3+i) 



(6.8) 



where A(j ;x,t; {yj}) is called the action. We shall as- 
sume that the force potential and the initial condition 
are periodic in the space variable and the period is taken 
to be unity. This assumption is very important for the 
discussion below. 

For a given initial condition at tj we next define a "min- 
imizing sequence" associated to (x, t) as a sequence of 
yi s (j = jo, jo + 1 , . . . , J(t)) at which the right-hand 
side of (6.7) achieves its minimum. Differentiating the 
action (6.8) with respect to the yj's one gets necessary 
conditions for such a sequence, which can be written as 
a sequence of (Euler-Lagrange) maps 



v 3+i = v j+fj(Vj), 
Vj+i = y3+ v 3+i( t 3+i - W 
= yj + (vj + fj{yj)){tj+i 



where 



Vj 

Vi = " 



Vj-i 



ti - t 



3-1 



(6.9) 
(6.10) 

(6.11) 



These equations must be supplemented by the initial 
and final conditions: 



where tl>o(x) is the initial potential. 

Using the variational principle we obtain the following 
"minimum representation" for the potential in the limit 
of vanishing viscosity which relates the solutions at any 



Vj =u (y jo ), (6.12) 
x = yj + v J+ i(t-tj). (6.13) 

It is easily seen that u(x, t) = vj + \ = (x — yj)/(t — tj). 
Observe that the "particle velocity" Vj is the velocity of 
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the fluid particle which arrives at yj at time tj and which, 
of course, has remained unchanged since the last kick (in 
Lagrangian coordinates). Equation (6.9) just expresses 
that the particle velocity changes by fj(yj) at the the 
kicking time tj . 

Note that (6.9)-(6.10) define an area-preserving and (ex- 
plicitly) invertible map. 

As in the case of continuous-in-timc forcing we can for- 
mulate the Burgers equation in the half-infinite time in- 
terval (— oo,t] without fully specifying the initial con- 
dition Uo(x) but only its (spatial) mean value (u) = 

Jq 1 u (x)dx. 

The construction of the solution in a half-infinite time 
interval is done by extending the concept of minimizing 
sequence to the case of dynamics starting at to = — oo. 
For a half- infinite sequence {yj} (j < J), let us define the 
action A({yj}; x, t; — oo) by (6.8) with j = — oo. Such 
a half-infinite sequence will be called a "minimizcr" (or 
"one-sided minimizer" ) if it minimizes this action with 
respect to any modification of a finite number of yj's. 
Specifically, for any other sequence {yj} which coincides 
with {yj} except for finitely many j's (i.e. yj — yj, j < 
J — k, k > 0), we require 

A({yj}; x,t;J-k)> A({y 3 }; x,t;J-k). (6.14) 



Of course, the Euler-Lagrange relations (6.9)-(6.10) still 
apply to such minimizcrs. Hence, if for a given x and t we 
know u(x, i), we can recursively construct the minimizer 
{yj} backwards in time by using the inverse of (6.9)- 
(6.10) for all j < J and the final condition - now an 
initial condition - (6.13) with vj + i — u(x,t). This is 
well defined except where u(x,t) has a shock and thus 
more than one value. 

One way to construct minimizers is to take a sequence 
of initial conditions at different times to — > — oo. At each 
such time some initial condition Uo(x) is given with the 
only constraint that it have the same prescribed value 
for (u). Then, (finite) minimizing sequences extending 
from t to t are constructed for these different initial 
conditions. This sequence of minimizing sequences has 
limiting points (sequences themselves) which are pre- 
cisely minimizers (E et al. 1998). The uniqueness of such 
minimizers, which would then imply the uniqueness of 
a solution to the Burgers equation in the time interval 
] — oo, t], can only be shown by using additional assump- 
tions, for example for the case of random forcing or when 
the forcing is time-periodic. 

If (u) = 0, the sequence {yj} minimizes the action 
A({yj}; x, t; — oo) in a stronger sense. Consider any 
sequence {yj} such that, for some integer P we have 
Vj = Vj + ^> J — J — k,k > and which differs ar- 
bitrarily from {yj} for j > J — k. (In other words, in 



a sufficiently remote past the hatted sequence is just 
shifted by some integer multiple of the spatial period.) 
We then have 



A({yj};x,t; -oo) > A({yj},x,t; -oo). 



(6.15) 



Indeed, for (u) = 0, the velocity potential for any initial 
condition is itself periodic. In this case a particle can be 
considered as moving on the circle S 1 and its trajectory 
is a curve on the space-time cylinder. The j/j's are now 
defined modulo 1 and can be coded on a representative 
< yj < 1. The Euler-Lagrange map (6.9)-(6.10) is still 
valid provided (6.10) is defined modulo 1. 

The condition of minimality implies now that yj and 
yj+i are connected by the shortest possible straight seg- 
ment. It follows that |uj+i| = p(yj,yj+i)/(tj + i — tj), 
where p is the distance on the circle between the points 
yj, yj+i, namely p(a, b) = min{|a — b\, 1— \a— b\}. Hence, 
the action A can be rewritten in terms of cyclic variables: 



A({y 3 };x,t;-co) 

+E 

3<J 



_ P 2 (x,yj) 
' 2(t-tj) 

2(tj+i - tj) 



F j+i(y. 



•3+1, 



. (6.16) 



The concept of "global minimizers" can be defined in 
a usual way. Namely, global minimizers correspond to 
one-sided minimizers that can be continued to a bilat- 
eral sequence {yj,— 00 < j < +00} while keeping the 
minimizing property. Such global minimizers correspond 
to trajectories of fluid particles that, from t = —00 to 
t = +00, have never been absorbed in a shock. As before 
we define a "main shock" as a shock which has always 
existed in the past. 

From now on we shall consider exclusively the case where 
the kicking is periodic in both space and time. Specifi- 
cally, we assume that the force in the Burgers equation 
is given by 



f(x,t)=g(x) J2 S (t-JT), 



g(x) = -^-G{x), 

ax 



(6.17) 
(6.18) 



where G(x), the kicking potential, is a deterministic 
function of x which is periodic and sufficiently smooth 
(e.g. analytic) and where T is the kicking period. The 
initial potential ip'mit (x) is also assumed smooth and pe- 
riodic. This implies that the initial velocity integrates to 
zero over the period. The case where this assumption is 
relaxed will be considered later in connection with the 
Aubry-Mather theory. 

The numerical experiments of [9] reported here have 
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Main Shock 



x 2n 

Fig. 28. Snapshots of the velocity for the unique time-peri- 
odic solution corresponding to the kicking force g(x) shown 
in the upper inset; the various graphs correspond to six out- 
put times equally spaced during one period. The origin of 
time is taken at a kick. Notice that during each period, two 
new shocks are born and two mergers occur. (From [9].) 

been made with the kicking potential 

G(x) — — sin3x + cos a;, (6.19) 

and a kicking period T= 1. Other experiments were done 
with (i) G(x) = — cosx and (ii) G(x) — (1/2) cos(2x) — 
cosx. The former potential produces a single shock and 
no preshock. As a consequence it displays no —7/2 law 
in the PDF of gradients. The latter potential gives es- 
sentially the same results as reported hereafter but has 
an additional symmetry. To avoid non-generic behaviors 
that could result from this symmetry it was chosen to 
focus on the forcing potential given by (6.19). 

The number of collocation points chosen for such simu- 
lations was mostly N x — 2 17 w 1.31 x 10 5 , with a few 
simulations done at N x = 2 20 (for the study of the re- 
laxation to the periodic regime presented below). Since 
the numerical method allows going directly to the de- 
sired output time (from the nearest kicking time) there 
is no need to specify a numerical time step. However, in 
order to perform temporal averages, e.g. when calculat- 
ing PDF's or structure functions, without missing the 
most relevant events (which can be sharply localized in 
time) sufficiently frequent temporal sampling is needed. 
The total number of output times N t w 1000, is thus 
chosen such that the increment between successive out- 
put times is roughly the two-thirds power of the mesh 
(this is related to the cubic structure of preshocks, see 
section 2.3). 

Figure 28 shows snapshots of the time-periodic solution 
at various instants. It is seen that shocks are always 
present (at least two) and that at each period two new 
shocks are born at t*\ w 0.39 and t* 2 ~ 0.67. There is one 
main shock which remains near x = ir and which collides 
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Fig. 29. Evolution of shock positions during one period. The 
beginnings of lines correspond to births of shocks (preshocks) 
at times t*i and t*2', shock mergers take place at times t c i 
and t C 2- The "main shock", which survives for all time, is 
shown with a thicker line. 
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Fig. 30. Exponential relaxation to a time-periodic solution 
for three different initial velocity data as labeled. The hori- 
zontal axis gives the time elapsed since t = 0. (From [9].) 

with the newborn shocks at t c \ s=s 0.44 and t C 2 ~ 0.86. 
Figure 29 shows the evolution of the positions of shocks 
during one period. 

It was found that, for all initial conditions uq(x) used, the 
solution u(x, t) relaxes exponentially in time to a unique 
function Uoo (x, t) of period 1 in time. Figure 30 shows the 

variation of J Q 71 \u(x, n_) — u ryo (x, 1_)| dx/(2n) for three 
different initial conditions as a function of the discrete 
time n. 

The phenomenon of exponential convergence to a unique 
space- and time-periodic solution is something quite gen- 
eral: whenever the kicking potential G(x) is periodic and 
analytic and the initial velocity potential is periodic (so 
that the mean velocity (u) ~0 at all times), there is ex- 
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poncntial convergence to a unique piecewise analytic so- 
lution. This can be proved rigorously (see Appendix to 
[9]) in the case when the functions G(x) have a unique 
point of maximum with a non-vanishing second deriva- 
tive (Morse generic functions). Here, we just explain the 
main ideas of the proof and give some additional prop- 
erties of the unique solution. 

One very elementary property of solutions is that, for 
any initial condition of zero mean value, the solution 
after at least one kick satisfies 

\u(x,t)\ < (1/2) + max|dG(a;)/dx|. (6.20) 

X 

Indeed, at a time t = n_ just before any kick we have 
x = y+u(x, n_) where y is the position just after the pre- 
vious kick of the fluid particle which goes to x at time n_ . 
It follows from the spatial periodicity of the velocity po- 
tential that the location y which minimizes the action is 
within less than half a period from x. Thus, \u(x, < 
1/2. The additional max x \dG(x)/dx\ term comes from 
the maximum change in velocity from one kick. Hence 
the solution is bounded. Note that if the spatial and tem- 
poral periods are L and T, respectively, the bound on 
the velocity becomes Lj (2T) + maxj |dG(x)/dx|. 

The convergence at large times to a unique solution can 
be understood in terms of the two-dimensional conserva- 
tive (area-preserving) dynamical system defined by the 
Euler-Lagrange map (6.9)-(6.10). By construction, we 
have u(x, 1+) = u(x) — dG(x)/dx, where u(x) is the so- 
lution of the unforced Burgers equation at time t = 1_ 
from the initial condition u(x) at time t — 0+. The map 
u i ► u(x) + g(x), where g(x) = — dG{x) / dx, here de- 
noted B g , solves the kicked Burgers equation over a time 
interval one. The problem is to show that the iterates 
BgU converge as n — > oo to a unique solution. 

If it were not for the shocks it would suffice to consider 
the two-dimensional Euler-Lagrange map. Note that, 
for the case of periodic kicking, this map has an obvi- 
ous fixed point P, namely (x = x c ,v = 0), where x c is 
the unique point maximizing the kicking potential. It is 
easily checked that this fixed point is an unstable (hy- 
perbolic) saddle point of the Euler-Lagrange map with 
two eigenvalues A = 1 + c + \J c? + 2c and 1/A, where 
c = -dl x G{x c )/2. 

Like for any two-dimensional map with a hyperbolic 
fixed point, there are two curves globally invariant by 
the map which intersect at the fixed point. The first is 
the stable manifold T^ s \ i.e. the set of points which con- 
verge to the fixed point under indefinite iteration of the 
map; the second is the unstable manifold T^ u \ i.e. the 
set of points which converge to the fixed point under in- 
definite iteration of the inverse map, as illustrated in fig- 
ure 31(a). Any curve which intersects the stable manifold 
transversally (at the intersection point, the two curves 
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Fig. 31. (a) Sketch of a hyperbolic fixed point P with stable 
(r^) and unstable (r^) manifolds. The dashed line gives 
the orbit of successive iterates of a point near the stable 
manifold, (b) Unstable manifold r' u ' on the (x, «)-cylinder 
(the ^-coordinate is defined modulo 1) which passes through 
the fixed point P — (x c ,0). The bold line is the graph of 
u ao (x, 1_). The main shock is located at xi — x r . Another 
shock at xi corresponds to a local zig-zag of F^> between A 
and B. 

are not tangent to each other) will, after repeated appli- 
cations of the map, be pushed exponentially against the 
unstable manifold at a rate determined by the eigenvalue 
1/A. In the language of Burgers dynamics, the curve in 
the (x, v) plane defined by an initial condition uq(x) will 
be mapped after time n into a curve very close to the 
unstable manifold. In fact, for the case studied numeri- 
cally, 1/A ~ 0.18 is within one percent of the value mea- 
sured from the exponential part of the graph shown in 
figure 30. Note that if the initial condition uo(x) contains 
the fixed point, the convergence rate becomes (1/A) 
(even higher powers of 1/ A are possible if the initial con- 
dition is tangent to the unstable manifold) . 

The fixed point P is actually a very simple global min- 
imizer: (yj — x Cl Vj — 0) for all positive and negative 
j's. It follows indeed by inspection of (6.16) that any 
deviation from this minimizer can only increase the ac- 
tion; actually, this trajectory minimizes both the kinetic 
and the potential part of the action. Note that the cor- 
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Fig. 32. Minimizers (trajectories of fluid particles) on the 
(x, t)-cylinder. Time starts at — oo. Shock locations at t = Cl- 
are characterized by the presence of two minimizers (an in- 
stance is at xi). The main shock is at xi — x r . The fat line 
x — x c is the global minimizer. 

responding fluid particle is at rest forever and will never 
be captured by a shock (it is actually the only particle 
with this property). It is easy to see that any minimizer 
is attracted exponentially to such a global minimizer as 
t — ► — oo. Thus, any point (yj,Vj) on a minimizer belongs 
to the unstable manifold and, hence, any regular 
part of the graph of the limiting solution Uoo (x) belongs 
to the unstable manifold r( u ). This unstable manifold is 
analytic but can be quite complicated. It can have sev- 
eral branches for a given x (see figure 31(b)) and does 
not by itself define a single- valued function u 00 (x). The 
solution has shocks and is only picccwisc analytic. Con- 
sideration of the minimizers is required to find the po- 
sition of the shocks in the limiting solution: two points 
with the same x corresponding to a shock, such as A and 
B on figure 31(b) should have the same action. 

Finally, we give the geometric construction of the main 
shock, the only shock which exists for an infinite time. 
Since the eigenvalue A is positive, locally, minimizers 
which start to the right of x c approach the global min- 
imizer from the right, and those which start to the left 
approach it from the left. Take the rightmost and left- 
most points x r and xi on the periodicity interval such 
that the corresponding minimizers approach the global 
minimizer from the right and left respectively (see fig- 
ure 32). These points are actually identical since there 
cannot be any gap between them that would have min- 
imizers approaching the global minimizer neither from 
the right nor the left. The solution u OQ (x) has then its 
main shock at xi = x r . 

6.2 Connections with Aubry-Mather theory 

In the previous subsection, the study of the solutions 
to the periodically kicked Burgers equation was limited 



to initial conditions with a vanishing spatial average 
b. With a non-vanishing mean velocity b, which in the 
forced case cannot be eliminated by a Galilean invari- 
ance, many of the properties of the solutions described 
above are still valid. However the action now depends on 

b. Global minimizers {Vj S \ j € Z} exist in this case as 
well. However generically they are not unique and do not 
correspond to fixed points of the Euler-Lagrange map 
(6.9)-(6.10). A global minimizer now minimizes the ac- 
tion 



Aoo({yk}) = A{{y k }; +oo; -oo) 



k— — oo 



2T 



(Vk+i-Vk-b) -G(y k+1 ) 



(6.21) 



This action is exactly the potential energy associated 
to an infinite chain of atoms linked by elastic springs 
and embedded in a periodic potential, problem known 
as the Frenkel-Kontorova model [52] . The parameter b 
represents the equilibrium length I of the springs and 
the spatial period L of the external potential (see fig- 
ure 33) is equal to 1. A global minimizer of (6.21) rep- 




Fig. 33. Sketch of the Frenkel-Kontorova model for the equi- 
librium states of an atom chain in a periodic potential. 

resents an equilibrium configuration of this system. The 
properties of this equilibrium, or ground states are de- 
termined by the competition between two tendencies: on 
the one hand the atoms tend to stabilize at those loca- 
tions where the potential is minimum; on the other hand, 
the springs tend to maintain them at a fixed distance of 
each other. When 6 = this competition disappears and 
the equilibrium is given by y k = x c , where x c is the loca- 
tion at which G attains its global minimum. For b ^ 0, 
the situation is more delicate and the structure of the 
ground states involves, as we shall now see, a problem 
of commensurate- incommensurate transition. The prop- 
erties of ground states were studied in great details by 
Aubry [3] and Mather [86]. 

The relations between the Burgers equation with a 
time-periodic forcing and Aubry-Mather theory were 
discussed for the first time in [70] and in [38] . The the- 
ory was further developed in [36,106]. For integer values 
of 6, the global minimizer is trivially associated to the 
fixed point (x,v) = (x Cl b) of the Euler-Lagrange map 
(6.9)-(6.10), which corresponds to a fluid trajectory lo- 
cated at integer times at x — x c and which moves on 
distance of b spatial periods during one temporal period. 
A similar argument implies that it is enough to study 
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values of b in the interval [0, l).To each global minimizcr 
{Vj B ^ j J ^ Z} is associated a rotation number defined as 

3=0 

which represents the time-average velocity of the mini- 
mizer. For a fixed value of the spatial average b of the 
velocity, all global minimizers associated to the solution 
of the Burgers equation have the same rotation number 
p. Indeed, as the dynamics is restricted to a compact do- 
main of the configuration space (in our case T), two min- 
imizers with different rotation numbers necessarily cross 
each other; this is an obvious obstruction to the action 
minimization property. In the case of rational rotation 
numbers the global minimizers correspond to periodic 
orbits of the dynamical system defined by the Euler- 
Lagrange map. An important feature is that for rational 
p, the rotation number does not change when varying 
b over a certain closed interval [6 m i n ,6 max ], called the 
mode-locking interval. On the contrary, irrational p cor- 
respond to a unique value of the parameter b. Such "ir- 
rational" values of b form a Cantor set of zero Lebesgue 
measure. In particular, the graph of p as a function of 
the parameter & is a "Devil staircase" (see figure 34) . 
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Fig. 34. Rotation number p as a function of the spatial mean 
of the velocity b for the standard map. 

When p is rational (p = p/q in irreducible form), global 
minimizers correspond to a periodic orbit of period q. 
It is easy to see that such an orbit generates q different 
but closely related global minimizers. Of course each of 
these global minimizer is the image of another one by the 
Euler-Lagrange map and is mapped back to itself after 
q iterations. This procedure generates a periodic orbit, 
which turns out to be hyperbolic one. Hence, each of 
the q global minimizers has a one-dimensional unstable 
manifold associated to it. The solution to the Burgers 
equation is formed by branches of these various mani- 
folds with jumps between them defining q global shocks. 



The picture is very different for values of b corresponding 
to irrational rotation numbers. Consider velocities and 
positions of all global minimizers at a fixed moment of 
time, say t = 0. They form a subset Q of the phase space 
TxB. Then two cases have to be distinguished: 

• The set Q forms a closed invariant curve for the Euler- 
Lagrange map. This invariant curve has a one-to-one 
projection onto the base T and dynamics on the curve 
is conjugated to a rigid rotation by angle p. The lim- 
iting solution of the Burgers equation is given by the 
invariant curve and does not contain any shocks. 

• The set Q forms an invariant Cantor set and the limit- 
ing solution of the Burgers equation contains an infi- 
nite number of shocks, none of which is a main shock. 

The Kolmogorov [80], Arnold [2] and Moser [93] the- 
ory (frequently referred to as KAM) describes invari- 
ant curves (or invariant tori) for small analytic per- 
turbations of intcgrable Hamiltonian systems, and thus 
the various types of dynamical trajectories. The KAM 
theory ensures that for sufficiently small perturbations, 
most of the invariant curves associated to Diophantine 
irrational rotation numbers are stable with respect to 
small analytic perturbations of the system. Diophantine 
irrational numbers possess fast converging approxima- 
tions by rational numbers (in a suitable technical sense). 
However, these invariant curves may disappear from the 
perturbed system when an interaction corresponding to 
a non-integrable perturbation gets sufficiently strong. 
Aubry Mather theory provides another variational de- 
scription for the KAM invariant curves. But even more 
importantly, it describes the invariant Cantor sets that 
appear instead of invariant curves in the case of strong 
nonlinear interactions. We have mentioned already that 
these Cantor sets correspond to global minimizers. Thus 
Aubry-Mathcr theory provides information about the 
global minimizers and, hence, allows one to study in such 
a situation the properties of limiting entropic solutions 
and, in particular, the structure of shocks. 

A numerical study of the Burgers equation in the inviscid 
limit, with periodic forcing and a non-vanishing spatial 
average of the velocity, reveals the appearance of shock 
accumulations. Such events occur for the values of the 
mean velocity b near the end-points of the mode-locking 
intervals, corresponding to rational rotation numbers. 
The shock accumulation phenomenon is due to the fact 
that the end-points 6 m in,&max of the mode- locking in- 
tervals can be approximated by convergent sequences of 
"irrational" values of the parameter b. This implies ac- 
cumulation of shocks, since for irrational rotation num- 
bers the number of shocks is infinite. 

The limiting solution u^ix, t) is completely determined 
by the function u{x) defined in the previous subsection. 
The function u{x) corresponds to a stroboscopic section 
of right after each impulse. The regular parts of u 
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are made of single- valued functions related to the unsta- 
ble manifolds. The shocks correspond to jumps, either 
between different branches of the same manifold (sec- 
ondary shocks) , or between the manifolds associated to 
different global minimizcrs (main shocks). 

When the rotation number is rational (p = p/q), there 
are q global minimizers. The positions of the q main 
shocks of u are determined by a requirement that the 
area defined by the graph of the solution is equal to the 
conserved quantity b. The latter constraint shows that 
the values of b compatible with the rotation number p/q 
belong to an interval [6 m i n , 6 ma x] bounded by the mini- 
mum and maximum areas defined by the unstable man- 
ifolds, as illustrated in figure 35. The detailed shape of 




Fig. 35. Sketch of the unstable manifolds of the two global 
minimizers associated to the rotation number p — 1/2. The 
values fe m i n and 6 max given by this configurations are repre- 
sented as grey areas. 

the manifolds can actually not be sketched on a figure. 
Generically the unstable manifold of a global minimizcr 
corresponding to a particular point of the basic periodic 
orbit of period q intersects transversally with the stable 
manifold of another minimizer corresponding to another 
point of the periodic orbit. Such an intersection leads to 
formation of a heteroclinic tangle, a notion which can be 
traced back to the work of Poincare. The heteroclinic in- 
tersection results in the formation of an infinite number 
of zig-zags of the unstable manifolds. These zig-zags are 
accumulating along the stable manifold and come arbi- 
trary close to the corresponding point of the periodic or- 
bit. The zig-zags contract exponentially in one direction 
(along the stable manifold) and are stretched exponen- 
tially in the other direction. It is easy to see that the 
accumulation of zig-zags generates an infinite number of 
"potential" shocks of smaller and smaller size which also 
accumulate near the periodic orbit. When the param- 
eter b is located well inside the mode-locking interval, 
the position of the main shock cuts off the accumulated 
shocks of small size so that the total number of shocks 
is of the order of unity. However, when b gets closer and 
closer to 6 max or 6 m i n , the main shocks move closer to 
the periodic points and a larger number of the small ac- 
cumulating shocks appears in the solution. This mecha- 
nism leads to an infinite number of shocks in the solution 
when b is equal to b min or 6 max (see figure 36(a)). Both 
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Fig. 36. (a) Accumulations of shocks occurring for b — 6 min or 
b = ftmax, due to the presence of an infinite number of loops 
of the unstable manifold in the homocline or heterocline 
tangle, (b) Shock accumulation at the fixed point (0, 0) of 
the standard map. Here, A = 0.1 and b = 0.15915. The latter 
value is close to the upper bound of the interval associated 
to p = 0. The upper inset is a zoom near (0, 0), illustrating 
the accumulation of shocks. 

the distances between two consecutive shocks and the 
sizes of the shocks decrease exponentially fast with the 
number of shocks; the rate is given by the stable eigen- 
value associated to the hyperbolic periodic orbit. It is 
interesting to mention that when b — b rn i n or b = b max 
the main shocks merge with the periodic orbit associ- 
ated to the global minimizers. Hence, for the end-points 
of the mode-locking interval the main shocks disappear. 

To illustrate numerically the change in behavior of the 
solution to the Burgers equation when the mean velocity 
b changes, we focus here on the simple periodic kicking 
potential G(x) — (X/2tt) cos(27rx) where A is a free pa- 
rameter. The associated Euler-Lagrange map then reads 

T: (y,v)\-^(y+v+\sm(2iry),v+\sm(2ny)). (6.23) 

This transformation is usually called the standard map 
(or Chirikov-Taylor map) . It is one of the simplest model 
for studying the presence of chaos in Hamiltonian dy- 
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Fig. 37. General aspect in position-velocity phase space of 
the dynamical system defined by the standard map (6.23) 
for two different values of the parameter (a) A = 0.1 and 
(b) A = 0.3. The corresponding time-periodic solutions to 
the kicked Burgers equation are represented as bold lines in 
both cases. The results are presented for the spatial mean 
velocities b — 0, b — 0.3 and b = 0.5. 

namical systems and in particular particularly to study 
the KAM theory. 

Figure 36(b) illustrates the accumulation of shocks due 
to the homoclinic or heteroclinic tangling for the first 
transition (starting from b = 0). This transition cor- 
responds to a rotation number of the global minimizcr 
changing value from p = to p > 0. When b is increased 
and gets close to the critical value, shocks accumulate 
on the left-hand side of the global minimizer located at 
(y,v) = (0,0). 

Other numerical experiments were performed in order 
to observe the destruction of invariant curves and the 
accumulation of shocks on Cantor sets for irrational ro- 
tation numbers, ft is of course impossible numerically to 
set the rotation number to an irrational value. Indeed, 
the values of b for which p is irrational are in a Cantor 



set. It is however possible to be very close to irrational 
rotation numbers. Figure 37 illustrates the changes in 
the behavior of the solutions to the periodically kicked 
Burgers equation when varying the parameter A. The 
time-asymptotic solutions associated to various values 
of the mean velocity b are shown for A = 0.1 and A = 0.3. 
For the latter value, all KAM invariant curves have al- 
ready disappeared. For 6 — and for all values of A the 
global minimizer trivially corresponds to the fixed point 
(0, 0) with a vanishing rotation number. For b = 0.5 
there are two global minimizers associated to the ratio- 
nal rotation number p — 1/2. For A = 0.1 and b = 0.3 
the rotation number is much closer to an irrational than 
in previous cases. The solution is then very close to the 
invariant curve associated to this value. Note that the 
main shock is actually located close to x w 0.85. It is so 
small that it can hardly be seen. When A = 0.3 the value 
b = 0.3 of the mean velocity no more corresponds to a 
rotation number close to an irrational value; it is now 
in the mode- locking interval associated to p = 1/3. This 
change in the rotation number reflects the dependence of 
the mode- locking intervals [6 m i n , b maK ] on the parameter 
A. The interval of values of b associated to p = is rep- 
resented as a function of A in figure 38. Such a structure 
is frequently called an Arnold tongue (see, e.g., [72]). 
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Fig. 38. Evolution as a function of the parameter A of the 
mode-locking interval [fe m i n ,6 m ax] associated to the rotation 
number p — 0. Such a graph is frequently referred to as an 
Arnold tongue. 

Finally, we discuss the structure of shocks in the case 
when the global minimizers form a Cantor set. There are 
then infinitely many gaps with no global minimizers. It 
is known in this case that all the gaps can be split into 
the finite number of images of the main gaps. For the 
standard map there is only one main gap. Its end-points 
(x 1 , v 1 ) and (x 2 , v 2 ) belong to the Cantor set associated 
to the global minimizers. All other gaps can be obtained 
by iterating this main gap with the Euler-Lagrange map 
(Standard map) for both positive and negative times: 
{xj,vj) = T^xV 1 ), (xlvf) = T J (x 2 ,i; 2 ), i e Z. One 
can show that the length of the ith gap tends to zero as 
i — ► ±oo. Since global minimizers are hyperbolic trajec- 
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tories one can connect the end-points of the main gap 
by two smooth curves: the stable manifold r( s ' and the 
unstable manifold r( u ) . As i — > oo the iterates of the sta- 
ble manifold T l T^ tend to a straight segment connect- 
ing the i-th gap with end-points at (xj, v}) and (x?,vf). 
The same is true for iterates of the unstable manifold 
7"«p( u ) i n the limit i — > — oo. On the contrary, negative 
iterates of the stable manifold and positive of the unsta- 
ble one form exponentially long curves connecting corre- 
sponding gaps. As usual we are interested in the iterates 
of the unstable manifold since they appear in the time- 
periodic solution of the Burgers equation. Such a solu- 
tion is formed by the iterates of the unstable manifold 
connecting all the gaps. Note that in the case of large 
negative i, the unstable manifold is close to a straight 
segment; hence there are no shocks located inside the 
corresponding gaps. Conversely, for large positive i, the 
unstable manifold is exponentially long and possesses 
large zig-zags. Hence, the solution to the Burgers equa- 
tion has one or several shocks inside such gaps. Since 
there are no shocks for gaps with large enough negative 
i, it follows that all the shocks have a finite age. In other 
words, the time-periodic solution has no main shocks. 

At the moment it was not possible to study numeri- 
cally the strange behavior of the solutions to the Burgers 
equation corresponding to global minimizers living on 
Cantor-like sets. Looking for such cases requires a very 
high spatial resolution in order to minimize the numeri- 
cal error in the approximation of the solution. Moreover, 
a large number of values for the parameters b and A has to 
be investigated in order to observe such a phenomenon. 
This would require heavy computer ressources. How- 
ever, many other aspects of the Aubry-Mather theory 
for Hamiltonian systems can be studied numerically us- 
ing the Burgers equation with periodic kicks. For in- 
stance it could be very useful for analyzing the higher 
dimensional versions. 

7 Velocity statistics in randomly forced Burgers 
turbulence 

The universality of small-scale properties in fully devel- 
oped Navier-Stokes turbulence has frequently been in- 
vestigated, assuming that a steady state is maintained 
by an external large-scale forcing. It is generally conjec- 
tured that the velocity increments have universal sta- 
tistical properties with respect to such a force. Under- 
standing this issue in simpler models of turbulence has 
motivated much work for over ten years. A toy model 
which has been extensively studied is the passive trans- 
port of a scalar field by random flows (see, e.g., [46]). 
Tools borrowed from statistical physics and field theory 
were used to describe and explain the anomalous scaling 
laws observed in the scalar spatial distribution. It was 
shown that the scale invariance symmetry is broken by 
geometrical constraints on tracer configurations that are 
statistically conserved by the dynamics. Universality of 



the intermittent scaling exponents with respect to the 
forcing was proven for the case where energy is injected 
at large scales [31,57,103,14]. 

Issues of universality for the nonlinear Burgers turbu- 
lence model has also been very much on the focus. The 
possibility to solve exactly a hydrodynamical problem 
displaying the same kind of quadratic nonlinearity as 
Navier-Stokes turbulence constitutes of course the cen- 
tral motivation. Three independent approaches were 
published almost simultaneously in 1995 and were at 
the origin of the growing interest in Burgers turbulence. 
First, an analogy was made in [22] between forced Burg- 
ers turbulence and the problem of a directed polymer in 
a random medium. This analogy was used to show that 
the shocks appearing in the solution lead to anomalous 
scaling laws for the structure functions. The strong 
intermittency could be related to the replica-symmetry- 
breaking nature of the disordered system associated 
to Burgers turbulence. This approach is discussed in 
subsection 7.1. Second, ideas using operator product 
expansions borrowed from quantum field theory were 
proposed in [99]. The goal was to close in the inertial 
range the equations governing the correlations of the ve- 
locity field in one dimension. This treatment of the dis- 
sipative anomaly is described in subsection 7.2. It yields 
a prediction for the probability density function (PDF) 
of velocity increments and gradients and in particular 
to a power-law behavior for the PDF of d x v at large 
negative values [99] . However, the value of the exponent 
of this algebraic tail has been a matter of controversy. 
An overview of the various works related to this issue is 
given in subsection 7.3. Finally, the turbulent model of 
the one-dimensional Burgers equation with a self-similar 
forcing was proposed in [30] as one of the simplest non- 
linear hydrodynamical problem displaying multiscaling 
of the velocity structure function. As stressed in subsec- 
tion 7.4 this problem is easily tractable numerically and 
some of the numerical observations can be confirmed by 
a one-loop renormalization group expansion. 

In what follows we consider the solutions to the Burg- 
ers equation with a homogeneous Gaussian random forc- 
ing that is delta-correlated in time. Namely, the spatio- 
temporal correlation of the forcing potential is taken to 
be 

(F(x,t)F(x',t')) = B{x-x')8{t-t'). (7.1) 

The function B contains information on the spatial 
structure of the forcing. It can be either smooth (i.e. 
concentrated at large spatial scales) or asymptoti- 
cally self-similar (i.e. behaving as a power law at small 
separations). In the former case the solution reaches 
exponentially fast a statistically stationary regime in 
any space dimension. The construction of the solution 
in this regime in terms of global minimizer and main 
shock is described in detail in section 5. When B does 
not decrease sufficiently fast at small separations (e.g. 
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B(r) ~ r 2h with h < 1 as r — > in one dimension), there 
is no rigorous proof of the existence of a statistically sta- 
tionary regime. However we assume in the sequel that 
such a stationary regime exists in order to perform a 
statistical analysis of the solutions to Burgers equation. 
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7.1 Shocks and bifractality - a replica variational ap- 
proach 

The replica solution for Burgers turbulence proposed 
in [22] is based on its analogy with the problem of a di- 
rected polymer in a random medium. As already stated 
in the Introduction, the viscous Burgers equation forced 
by the potential F is equivalent to finding the partition 
function Z of an elastic string in the quenched spatio- 
temporal disorder V(x, t) = F(x,t)/2v (remember that 
t has to be interpreted as the space direction in which 
the polymer is oriented) . This relation is obtained by ap- 
plying to the velocity potential ^ the Hopf-Cole trans- 
formation Z(x, t) = exp(ty(x, t)/2v). The solution of the 
problem can be written in terms of the path integral 



Z(x,t) = [ exp(-W(7)) d[7(-)] , 
l 



with H(j) 



2v 



%s) +F(7(s),s) 



ds. (7.2) 



In the analogy between Burgers turbulence and directed 
polymers, the polymer temperature is assumed to be 
unity and its elastic modulus is l/(2u). The strength 
of the potential fluctuations applied to the polymer de- 
pends on the viscosity and is oc e x l 2 Lf j(2v) (where e is 
the energy injection rate and Lf is the spatial scale of 
forcing) . In order to calculate the various moments of the 
velocity field v = — V^, one needs to average the loga- 
rithm of the partition function Z, a celebrated problem 
in disordered systems. 

Bouchaud, Mezard and Parisi proposed in [22] the use 
of a replica trick in order to estimate the average free 
energy (InZ). The first step is to write the zero- replica 
limit InZ = lim„^o {Z n — l)/n. Then, the moments 
(Z n ) are used to generate an effective attraction between 
replicas: they are written as the partition functions of 
the disorder- averaged Hamiltonian TL n {li , • ■ • > 7n) asso- 
ciated to n replicas of the same system [90] 
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where B denotes the spatial part of the forcing potential 
correlation. The next step is to study this problem by a 
variational approach. The Hamiltonian H n is replaced 
by an effective Gaussian quadratic Hamiltonian that can 



The kernel Qij is then chosen in such a way that it 
minimizes the free energy. It is shown in [22] that the 
optimal Gaussian Hamiltonian is the solution of a sys- 
tem of equations that can be solved following the ansatz 
proposed in [89]. When d > 2 this approach singles 
out two regimes depending on the Reynolds number 
Re = e 1 ^ 3 ^ 4 jv. These regimes are separated by the 

critical value Re c = [2(l-2/c0 1_<i/2 ] 1/3 - When Re < Re c 
the optimal solution is of the form Qij = Q Sij + Q\ and 
obeys the replica symmetry. In finite-size systems it cor- 
responds to a linear velocity profile. When Re > Re c 
the correct solution is given by the one-step replica- 
symmetry-breaking scheme (see [89]). The off-diagonal 
elements of Qij are then parameterized with two func- 
tions depending on whether the indices i and j belong to 
the same block or to different blocks. Qualitatively, the 
one-step replica-symmetry-breaking approach amounts 
to the assumption that the instantaneous velocity po- 
tential can be written as a weighted sum of Gaussians, 
leading to an approximation of the velocity field as 

2v Y< a (x-r a )e~ ReUciHS ~ ?A2,2L2f) 
v(x,t)~— a -ReU.+UZ-P.WWf) ' (7 ' 5) 

Z—ia 

where the f a 's are independent variables with a Poisson 
distribution of density exp(— /). The r a are uniformly 
and independently distributed in space. In (7.5) the sum 
over a is running from 1 to a large-enough integer M. 
The typical shape of the approximation of the veloc- 
ity field given by (7.5) is represented in figure 39(a) in 
the two-dimensional case. In the limit of large Reynolds 
numbers the random velocity field given by (7.5) typi- 
cally contains cells of width oc Lf. The width of a shock 
separating two cells is of the order of Lf / Re. 

The replica approximation (7.5) leads to an estimate of 
the PDF p(Av, r) of the longitudinal velocity increment 
Av = (v(x + r e, t) — v(x, t)) • e , where e is an arbitrary 
unitary vector. When Re » 1 and r <C Lf this approxi- 
mation takes the particularly simple asymptotic form 



p(Av,r) w S[ Av - U f 



r 1 /Av\ 



where Uf = RevjLf is the typical velocity associated 
to the scale L f and g is a scaling function that is deter- 
mined explicitly in [22] . This approximation is in agree- 
ment with the following qualitative picture. With a prob- 
ability almost equal to one, the two points x and x + re 
lie in the same cell; the velocity increment is then given 
by the typical velocity gradient which, according to the 
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Fig. 39. (a) Typical shape of the velocity field given by the 
replica approximation in dimension d = 2 obtained from 
(7.5) for Re — 10 3 . The contour lines represent the velocity 
modulus. Note the cell structure of the domain, (b) Scaling 
exponents of the pth order structure function. 

approximation (7.5), is order Uf/Lf. With a probabil- 
ity r/Lf the two points are sitting on different sides of 
a shock separating two such cells and the associated ve- 
locity difference is of the order of Uf. 

The structure functions of the velocity field given by 
the various moments of Av can be straightforwardly 
estimated from the approximation (7.6). Their scaling 
behavior (Av p ) ~ at small separations r display 
a bifractal behavior as sketched in figure 39(b). When 
p < 1, the first term on the right-hand side of (7.6) domi- 
nates and (Av p ) oc Uf (r/Lf) p . Forp > 1 the shock con- 
tribution is dominating the small-r behavior and thus 
(A^)oc U p f {r/Lf). 

This approach, which makes use of replica tricks, is as we 
have seen able to catch the leading scaling behavior of 
velocity structure functions in any dimension. It is based 
on approximations of the velocity field by the superposi- 
tion (7.5) of Gaussian velocity potentials. A first advan- 
tage of this method is that it catches the generic aspect 
of the solution including the hierarchy of high-order sin- 
gularities appearing in the solution when Re — > oo which 
was examined in section 2.3. This method also gives 
predictions regarding the dependence on Re of the sta- 



tistical properties of the solution. However, as stressed 
in [22], the validity of this approximation is expected to 
hold only in the limit of infinite space dimension d. In 
particular, it is known that for d < 2 a full continuous 
rcplica-symmetry-breaking scheme is needed [89] . Nev- 
ertheless, as we have seen, there is enough evidence that 
this approach describes very well the qualitative aspects 
of the solution. 

7.2 Dissipative anomaly and operator product expan- 
sion 

The replica-trick approach described in the previous 
subsection cannot reproduce one of the main statis- 
tical features of the solution, namely the tails of the 
velocity increments PDF p(Av,r). Indeed the predic- 
tion (7.6) based on a variational approximation of the 
velocity field implies that p identically vanishes when 
Av > Uf (r/Lf). In order to study the quantitative be- 
havior of the PDF p(Av,r) in the inviscid limit v — > 
(or cquivalcntly Re — ► oo), Polyakov [99] proposed to 
use an operator product expansion. This approach leads 
to an explicit expression for p(Av, r) and predicts a 
super-exponential tail at large positive values and a 
power-law behavior for negative ones. Such predictions 
have immediate implications for the asymptotics of the 
PDF p(£) of the velocity gradient £ = d x v. The work 
of Polyakov was the starting point of a controversy on 
the value of the exponent of the left tail of Before 
returning to this issue in the next subsection, we give 
in the sequel a quick overview of the original work by 
Polyakov. 

We henceforth focus on the one-dimensional solutions 
to the Burgers equation with Gaussian forcing whose 
autocorrelation is given by (7.1). Following [99] (see also 
[19,20]) we introduce the characteristic function of the 
n-point velocity distribution 

Z n {\j, Xj ;t) = (e*iv(xut)+-+\ n v(x n ,t)^ ( 77 ) 

For a finite value of the viscosity i^, it is easily seen that 
this quantity is a solution to a Fokker-Planck (master) 
equation obtained by differentiating Z n with respect to 
t and using the Burgers equation and the fact that the 
forcing is Gaussian and <5-correlated in time. This leads 
to 

dz n v d ( i dz n \ = 

dt If ■ ' W U, dxi) 

= ^ ^2 b(xi - Xj) Aj Xj Z n + , (7.8) 

where b = (d 2 B)/(dr 2 ) denotes the spatial part of the 
correlation of the forcing applied to the velocity field. 
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T>v %s> denotes the contribution of the dissipative term and 
reads 



2)( n ) 



E 



(7.9) 



This term does not vanish in the limit v — > since the 
solutions to the Burgers equation develop singularities 
with a finite dissipation. It has been proposed in [99] to 
use an analogy with the anomalies appearing in quantum 
field theory in order to tackle this term in the inviscid 
limit. The important assumption is then made that the 
singular term in the operator product expansion relates 
linearly to the characteristic function Z n . Since this ex- 
pansion should preserve the statistical symmetries of the 
Burgers equation, it leads to the replacement in all av- 
erages of the singular limit lim„^o v A (d 2 v) e Xv by the 
asymptotic expression 



a b- 1 d d 



(7.10) 



where the coefficients a, b and c are parameters that can 
be determined only indirectly. However their possible 
values can be restricted by requiring that Z n is the char- 
acteristic function of a probability distribution which is 
non-negative, finite, normalizablc, and that the dissipa- 

(n) 

tive term V v acts as a positive operator. Finding these 
coefficients is similar to an eigenvalue problem in quan- 
tum mechanics. 

We now come to a crucial point in Polyakov's approach. 
Important restrictions on the form of the different 
anomalous terms in (7.10) result from the fact that the 
solutions to the Burgers equation obey a certain form of 
Galilean invariance. A notion of "strong Galilean prin- 
ciple" is introduced for invariance of the n-point distri- 
bution of velocity under the transformation v i— » v + v 
with vo an arbitrary constant. As a consequence, the n- 
point characteristic function Z n has to be proportional 
to S(Xi + • • • + A„). The operators appearing in the limit 
v — > have to be consistent with such an invariance. 
In [99] it is argued that this symmetry is automatically 
broken by the forcing that introduces a typical velocity 
[v 2 ) 1 / 2 cx 6 1 / 3 (0)L 1 ' 3 . However Polyakov assumes this 
"strong Galilean principle" to be asymptotically recov- 
ered in the limit L — > oo of infinite-size systems. In the 
case of finite-size systems, when L is of the order of the 
correlation length Lf of the forcing, the strong Galilean 
symmetry is broken because of the conservation of the 
spatial average of v which introduces a characteristic 
velocity v n — (1/L) J v(x,t) dx. However, the Galilean 
symmetry should be recovered when averaging the cor- 
relation functions with respect to the mean velocity ^o- 
This symmetry restoration was introduced in [20] where 
it is referred to as the "weak Galilean principle". The 
n-point characteristic function associated to an aver- 



age velocity Vo relates to that associated to a vanishing 
mean velocity by 

Z n (Xj,Xj;t;v ) ^c" ^^ 3 Z n (\j, xj; t; 0) . 
After averaging with respect to u , one obtains 

Z n (Xj,Xj;t) = 2n6\ y^Aj Z n (Xj,Xj;t; 0) . (7.11) 



One can easily check that (7.8), together with the dis- 
sipative term given by (7.10), are compatible with this 
expression for the n-point characteristic function Z n . 
Moreover, any higher-order term in the expansion (7.10) 
of the dissipative anomaly would violate Galilean invari- 
ance. 

To obtain the statistical properties of the solution, one 
needs to further restrict the values of the three free pa- 
rameters a, b, and c appearing in the expansion (7.10). 
Following [99] this can be done by considering the case 
n = 2 that corresponds to the equation for the PDF of 
velocity differences. Performing the change of variables 
Ai ; 2 = A ± fj, and X\ t2 — X ± y/2, and assuming that 
A <C /i and y <C Lf (so that the spatial part of the 
forcing correlation is to leading order b(y) ~ &o — biy 2 ), 
the stationary and space-homogeneous solutions to the 
master equation (7.8)) satisfy 
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dfidy 



- (2b A 2 



bin 2 y 2 )Z 2 



2b dZ 2 , dZ 2 dZ 2 
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fi oy oA oil 



— -+cfi— -. (7.12) 



It is next assumed in [99] (see also [20]) that the velocity 
difference v{x\,t) — v(x 2 ,t) is statistically independent 
of the mean velocity (v(x\, t) + v(x 2 ,t))/2. This implies 
that the two-point characteristic function factorizes as 
Z 2 = Z 2 {A)Z 2 {fJ., y), where the two functions Z^ and 
Z^ satisfy the closed equations 



26 A 2 Z+ = cA^, 
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dfidy 



> -b^ 2 y 2 Z 2 = a Z 2 + 2bdZ2 



M dy 



dZ: 



(7.13) 
(7.14) 



The solution to the first equation corresponds to a Gaus- 
sian distribution which is normalizablc only if c < 0. As 
shown numerically in [20] this distribution is representa- 
tive of the bulk of the one-point velocity PDF. Informa- 
tion on the solutions to the second equation can be ob- 
tained assuming the scaling property Z% (u, y) — 
which amounts to considering only those contributions 
to the distribution of velocity differences stemming from 
velocity gradients £ = d x v. This yields a prediction the 
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negative and positive tails of the PDF of velocity gradi- 
ents: 

p(0 oc |CP Q when £-> -co , (7.15) 
p(0 cx C /3 exp(-Ce 3 ) when ^+00, (7.16) 

where C is a constant, which depends only on the 
strength of the forcing. The two exponents a and (3 are 
related to the coefficient b of the anomaly by 

a = 2b+land/? = 2b-l. (7.17) 

The value of b remains undetermined but is prescribed 
to belong to a certain range. This approach was first 
designed in [99] for infinite-size systems where strong 
Galilean invariance holds. In that case consistency with 
such an invariance leads to dropping the third term in 
the operator product expansion (i.e. c = 0). Positivity 
and normalizability of the two-point velocity PDF and 
non-positivity of the anomalous dissipation operator im- 
ply that the two other coefficients form a one-parameter 
family with 3/4 < b < 1. In particular, this implies 
that the left tail of the velocity gradient PDF with ex- 
ponent a should be shallower than £ -3 . As we will see in 
the next section, strong evidence has been obtained that 
p(£) ex £~ 7 / 2 for £ — > — oo. This seems to contradict the 
approach based on operator product expansion. How- 
ever, as argued in [20] , the breaking of strong Galilean in- 
variance occurring in finite-size systems and resulting in 
the presence of the c anomaly broadens the range of ad- 
missible values for b. In particular it allows for the value 
b = 5/4 which corresponds to the exponent a = 7/2. 

7.3 Tails of the velocity gradient PDF 

After the numerical work of Chekhlov and Yakhot [29] , 
the asymptotic behavior at large positive and negative 
values of the PDF of velocity derivatives £ = d x v for 
the one-dimensional randomly forced Burgers equation 
attracted much attention. A broad consensus emerged 
around the prediction of Polyakov [99] that p(£) dis- 
plays tails of the form (7.16) and (7.15), but the values 
of the exponents a and /3 were at the center of a contro- 
versy. Note that the presence of a super-exponential tail 
oc exp(— C£ 3 ) at large positive arguments has been con- 
firmed by the use of instanton techniques [60] and that 
the only remaining uncertainty concerns the exponent 
of the algebraic prefactor. A standard approach to de- 
termine the exponents a and f3 appearing in (7.15) and 
(7.16) makes use of the stationary solutions to the in- 
viscid limit of the Fokkcr Planck equation for the PDF, 
namely 

-ds(ep)-&+vds[(d 3 x v\d x v = Z)p] = b%p. (7.18) 

Here the brackets (-|-) denote conditional averages and 
the right-hand side expresses the diffusion of probability 



due to the delta-correlation in time of the forcing. The 
main difficulty in studying the solutions of (7.18) stems 
from the treatment of the dissipative term D v (£) = 
v ®£ [(9x v \9xV — £,) p] in the limit v ^> 0. Thevaluea = 3 
is obtained if a piecewise linear approximation is made 
for the solutions of the Burgers equation [21]. Gotoh and 
Kraichnan [59] argued that the dissipative term is to 
leading order negligible and presented analytical and nu- 
merical arguments in favor of a — 3 and (3=1. However, 
the inviscid limit of (7.18) contains anomalies due to the 
singular behavior of D l '(£) in the limit v — > 0. As we 
have seen in previous section, the approach based on the 
use of an operator product expansion [99] leads to a rela- 
tion involving unknown coefficients which must be deter- 
mined, e.g., from numerical simulations [111,19,20], and 
restricts the possible values to 5/2 < a < 3 [6]. Anoma- 
lies cannot be studied without a complete description of 
the singularities of the solutions, such as shocks, and a 
thorough understanding of their statistical properties. 

E, Khanin, Mazel and Sinai made a crucial observa- 
tion in [37] that large negative gradients stem mainly 
from preshocks, that is the cubic-root singularities in 
the velocity preceding the formation of shocks (see sec- 
tion 2.3). They then used a simple argument for de- 
termining the fraction of space-time where the veloc- 
ity gradient is less than some large negative value. This 
leads to a — 7/2, provided preshocks do not cluster. 
Later on, this approach has been refined by E and Van- 
den Eijnden who proposed to determine the dissipative 
anomaly of (7.18) using formal matched asymptotics [39] 
or bounded variation calculus [42]. As we shall see be- 
low, with the assumption that shocks are born with a 
zero amplitude, that their strengths add up during colli- 
sions, and that there ar no accumulations of preshocks, 
the value a = 7/2 was confirmed [42]. Other attempts to 
derive this value using also isolated preshocks have been 
made [81,6]. Note that there are simpler instances, in- 
cluding time-periodic forcing [9] (see section 6) and de- 
caying Burgers turbulence with smooth random initial 
conditions [8,42] (see section 4.1), which fall in the uni- 
versality class a — 7/2, as can be shown by systematic 
asymptotic expansions using a Lagrangian approach. 

We give here the flavor of the approach used in [39] 
in order to estimate the dissipative anomaly D°(£) = 

linwo-D^O- 0ne first notices that for |£| > 6 1 / 3 , the 
forcing term in the right-hand side of (7.18) becomes 
negligible, so that stationary solutions to the Fokker- 
Planck equation satisfy 

p(o«ier 3 / ^'WiO- aw) 

J — oo 

A straightforward consequence of this asymptotic ex- 
pression is that, if the integral in the right-hand side de- 
creases as £ — > — oo (i.e. if £D V (£) is integrable), then 
p(£) decreases faster than |£| -3 , and thus a > 3. 
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To get some insight into the behavior of D v as v —* 0, one 
next observes that the solutions to the one-dimensional 
Burgers equation contain smooth regions where viscos- 
ity is negligible, which are separated by thin shock layers 
where dissipation takes place. The basic idea consists in 
splitting the velocity field v into the sum of an outer so- 
lution away from shocks and of an inner solution near 
them for which boundary layer theory applies. Matched 
asymptotics are then used to construct a uniform ap- 
proximation of v. To construct the inner solution near 
a shock centered at x — x*, one performs the change of 
variable ihi=(i- x*)/v and looks for an expression 
ofv(x,t) = v{x+ + vx, t) in the form of a Taylor expan- 
sion in powers of v. v = vq + isvi +o(V). At leading order, 
the inner solution satisfies 

[vo-v*]d £ v Q =dlv Q , (7.20) 

where = (dx*)/(dt). This expression leads to the well- 
known hyperbolic tangent velocity profile 

s ( sx \ 
So =«*- -tanhl —J . (7.21) 

Here, s = v(x±+, t) — v(x+ — , t) denotes here the velocity 
jump across the shock and is given by matching condi- 
tions to the outer solution. The term of order v is then 
a solution of 

dm + [v - «*] a £ Si = d\v x + /(x, t). (7.22) 

In order to evaluate the dissipative anomaly, it is con- 
venient to assume spatial ergodicity so that the viscous 
term in (7.18) can be written as 

£>"(0 = ud t lrm -L j\x 8 3 x v 5(d x v-£). (7.23) 

L^oo IL J_ L 

In the limit v — > the only remaining contribution stems 
from shocks and is thus given by the inner solution. Using 
the expansion of the solution up to the first order in v. 
this leads to writing the dissipative term in the limit of 
vanishing viscosity as (see Appendix of [41] for details) 

D\i) = \ J\ss\p+(s,0+p-(s,0], (7-24) 

where p is the density of shocks and p + (respectively p~ ) 
is the joint probability of the shock jump and of the value 
of the velocity gradient at the right (respectively left) of 
the shock. This expression guarantees the fmiteness of 
the dissipative anomaly and in particular the fact that 
the integral in the right-hand side of (7.19) is finite in 
the limit v — > and converges to 0. As a consequence, 
this gives a proof that the exponent a of the left tail of 
the velocity gradient PDF is larger than 3. 



To proceed further, E and Vanden Eijnden proposed to 
estimate the probability densities p + and p~ by deriving 
master equations for the joint probability of the shock 
strength s, its velocity and the values ^ of the veloc- 
ity gradient at its left and at its right. This is done in [42] 
using a formulation of Burgers dynamics stemming from 
bounded variation calculus. More precisely, it is shown 
in [108] that the Burgers equation is equivalent to con- 
sidering the solutions to the partial differential equation 

8 t v + vd x v = f , (7.25) 

where v(x,t) = (v(x+,t) + v(x— , t))/2. Basically this 
means that Burgers dynamics can be formulated in 
terms of the transport of the velocity field by its average 
v. This formulation straightforwardly yields a master 
equation for v(x±, t) and d x v(x±, t) which is then used 
to estimate p^ and the dissipative anomaly (7.24). Al- 
though the treatment of the master equation does not 
involve any closure hypothesis, it is not fully rigorous: 
in particular it requires the assumption that shocks are 
created with zero amplitude and that shock amplitudes 
add up during collision. However such an approaches 
strongly suggests that a = 7/2 and j3 = 1. 

Obtaining numerically clean scaling for the PDF of gra- 
dients is not easy with standard schemes. Let us ob- 
serve that any method involving a small viscosity, either 
introduced explicitly (e.g. in a spectral calculation) or 
stemming from discretization (e.g. in a finite difference 
calculation), may lead to the presence of a power-law 
range with exponent —1 at very large negative gradi- 
ents [59]. This behavior makes the inviscid range 
appear shallower than it actually is, unless extremely 
high spatial resolution is used. In contrast, methods that 
directly capture the inviscid limit with the appropriate 
shock conditions, such as the fast Legendre transform 
method [94], lead to delicate interpolation problems. 
They have been overcome in the case of time-periodic 
forcing [9] but with white-noise-in-time forcing, it is dif- 
ficult to prevent spurious accumulations of preshocks 
leading to a = 3. 

To avoid such pitfalls, a Lagrangian particle and shock 
tracking method was developed in [6]. This method is 
able to separate shocks and smooth parts of the solution 
and is particularly effective for identifying preshocks. 
The main idea is to consider the evolution of a set of 
N massless point particles accelerated by a discrete-in- 
time approximation of the forcing with a uniform time 
step. When two of these particles intersect, they merge 
and create a new type of particle, a shock, characterized 
by its velocity (half sum of the right and left velocities 
of merging particles) and its amplitude. The particle- 
like shocks then evolve as ordinary particles, capture 
further intersecting particles and may merge with other 
shocks. In order not to run out of particles too quickly, 
the initial small region where particles have the least 
chance of being subsequently captured is determined by 
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localization of the global minimizcr of the Lagrangian 
action (see section 5.1). The calculation is then restarted 
from t = for the same realization of forcing but with 
a vastly increased number of particles in that region. 
This particle and shock tracking method gives complete 
control over shocks and preshocks. 
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Fig. 40. PDF of the velocity gradient at negative values in 
log-log coordinates obtained by averaging over 20 realiza- 
tions and a time interval of 5 units of time (after relaxation 
of transients) . The simulation involves up to N — 10 5 parti- 
cles and the forcing is applied at discrete times separated by 
St — 10~ 4 . Upper inset: local scaling exponent (from [6]). 

Figure 40 shows the PDF of the velocity gradients in 
log-log coordinates at negative values, for a Gaussian 
forcing restricted to the first three Fourier modes with 
equal variances such that the large-scale turnover time is 
order unity. Quantitative information about the value of 
the exponent is obtained by measuring the "local scaling 
exponent", i.e. the logarithmic derivative of the PDF 
calculated in this case using least-square fits on half- 
decades. It is seen that over about five decades, the local 
exponent is within less than 1% of the value a = 7/2 
predicted by E et al. [37]. 

7.4 Self-similar forcing and multiscaling 

As we have seen in section 7.1, the solutions to the Burg- 
ers equation in a finite domain and with a large-scale 
forcing have structure functions (moments of the ve- 
locity increment) displaying a bifractal scaling behav- 
ior. Such a property can be easily interpreted by the 
presence of a finite number of shocks with a size order 
unity in the finite system. Somehow this double scal- 
ing and its relationship with singularities gives some in- 
sight on the multiscaling properties that are expected 
in the case of turbulent incompressible hydrodynamics 
flows. There is a general consensus that the turbulent 
solutions to the Navier-Stokes equations display a full 
multifractal spectrum of singularities which are respon- 
sible for a nonlinear p-dependence of the scaling expo- 



nents Q p associated to the scaling behavior of the p-th 
order structure function [53] . The construction of simple 
tractable models which are able to reproduce such a be- 
havior has motivated much work during the last decades. 
Significant progress, both analytical and numerical, has 
been made in confirming multiscaling in passive-scalar 
and passive- vector problems (see, e.g., [46] for a review). 
However, the linearity of the passive-scalar and passive- 
vector equations is a crucial ingredient of these studies, 
so it is not clear how they can be generalized to fluid 
turbulence and the Navier-Stokes equation. 

After the work of Chckhlov and Yakhot [30] , it appeared 
that the Burgers equation with self-similar forcing could 
be the simplest nonlinear partial differential equation 
which has the potential to display multiscaling of veloc- 
ity structure functions. We report in this section various 
works that tried to confirm or to weaken this statement. 
Let us consider the solutions to the one-dimensional 
Burgers equation with a forcing term f(x, t) which is ran- 
dom, space-periodic, Gaussian and whose spatial Fourier 
transform has correlation 

(/>, t)f(k', t 1 )) = 2D \kf S(t - t') S(k + k'). (7.26) 

The exponent determines the scaling properties of the 
forcing. When > the force acts at small scales; for 
instance = 2 corresponds to thermal noise for the ve- 
locity potential, and thus to the KPZ model for interface 
growth [74]. It is well known in this case (see, e.g., [5]) 
that the solution displays simple scaling (usually known 
as KPZ scaling), such that £ g = q for all q. More gener- 
ally, the case > can be exactly solved using a one- 
loop renormalization group approach [88] . 

As stressed in [64] , renormalization group techniques fail 
when < and the forcing acts mostly at large scales 
and non-linear terms play a crucial role. When < — 3, 
the forcing is differentiable in the space variable, the so- 
lution is piecewise smooth and contains a finite number 
of shocks with sizes order unity. The scaling exponents 
are then £ p = min (1, p). In the case of non-diffcrcntiablc 
forcing (—3 < < 0), the presence of order-unity shocks 
and dimensional arguments suggest that the scaling ex- 
ponents are ( p = min (1, —p0/3). However, very little is 
known regarding the distribution of shocks with inter- 
mediate sizes. In particular, there is no clear evidence 
whether or not they form a self-similar structure at small 
scales. We summarize here some studies which were done 
on Burgers turbulence with self-similar forcing to show 
how difficult it might be to measure scaling laws of struc- 
ture functions and in particular how logarithmic correc- 
tions can masquerade anomalous scaling. 

For this we focus on the case = — 1 which has attracted 
much attention; indeed, dimensional analysis suggests 
that ( p = p/3 when p < 3, leading to a K41-type —5/3 
energy spectrum. Early studies [29,30] seemed to confirm 
this prediction using pseudo-spectral viscous numerical 
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Fig. 41. Representative snapshots of the velocity v (jagged 
line) in the statistically stationary regime, and of the inte- 
gral of the force / over a time step (rescaled for plotting 
purposes). 

simulations at rather low resolutions (around ten thou- 
sands gridpoints). It was moreover argued in [64,65] that 
a self-similar forcing with — 1 < (3 < 0, could lead to gen- 
uine multifractality. The lack of accuracy in the determi- 
nation of the scaling exponents left open the question of 
a weak anomalous deviation from the dimensional pre- 
diction. This question was recently revisited in [91] with 
high-resolution inviscid numerical simulations using the 
fast Legendre transform algorithm (see section 2.4.2). 
A typical snapshot of the forcing and of the solution in 
the stationary regime are represented in figure 41. It is 
clear that because of shocks the velocity develops small- 
scale fluctuations much stronger than those present in 
the force. However one notices that shock dynamics and 
spatial finitcncss of the system lead, as predicted, to the 
presence of few shocks with order-unity sizes. 




Fig. 42. Scaling exponents ( p versus order p for 
N = 2 16 (o), 2 18 (*), and 2 20 (o) grid points. Error bars (see 
text) are shown for the case N = 2 20 . The deviation of £ p 
from the exponents for bifractal scaling (full lines), shown 
as an inset, naively suggests multiscaling (from [91]) 

Structure functions were measured with high accuracy. 
They typically exhibit a power-law behavior over nearly 
three decades in length scale; this is more than two 
decades better than in [30]. In principle one expects to 
be able to measure the scaling exponents with enough 



accuracy to decide between bifractality and multiscal- 
ing. Surprisingly the naive analysis summarized in fig- 
ure 42 does suggest multiscaling: the exponents ( p de- 
viate significantly from the bifractal-scaling prediction 
(full lines). Since the goal here is to have a precise han- 
dle on the scaling properties of velocity increments, it is 
important to carefully define how the scaling exponents 
are measured. They are estimated from the average log- 
arithmic derivative of Sp hs (r) = (\v(x + r) — v(x) \ p ) over 
almost two decades in the separation r. The error bars 
shown are given by the maximum and minimum devia- 
tions from this mean value in the fitting range. Note also 
that the observed multiscaling is supported by the fact 
that there is no substantial change in the value of the 
exponents when changing the number N of grid points 
in the simulation from 2 16 to 2 20 : any dependence of 
C p upon N is much less than the error bars determined 
through the procedure described above. 
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Fig. 43. Log-log plots of Sl hs (r) (dashed line), S 3 (r) 
(crosses), and ((d + v)' i ) (squares) versus r. The continuous 
line is a least-square fit to the range of points limited by two 
vertical dashed lines in the plot. Inset: An explicit check of 
the von Karman-Howarth relation (7.27) from the simula- 
tions with N = 2 20 reported in [91]. The dashed curve is the 
integral of the spatial part of the forcing correlation and the 
circles represent the numerical computation of the left-hand 
side. 

As found in [91], the observed deviations of the scaling 
exponents from bifractality are actually due to the con- 
tamination by subleading terms in (r). To quantify 
this effect, let us focus on the third-order structure func- 
tion (p = 3) for which one measures (3 « 0.85±0.02 over 
nearly four decades (see figure 43). To estimate sublead- 
ing terms we first notice that the third-order structure 
function S^r) = ((v(x + r) — w(x)) 3 ), which is defined, 
this time, without the absolute value, obeys an analog of 
the von Karman-Howarth relation in fluid turbulence, 
namely 

±S 3 (r) = £b(r')dr', (7.27) 

where b(-) denotes the spatial part of the force corre- 
lation function, i.e. (f(x + r,t')f(x,t)) = b(r)S(t — t'). 
This relation, together with the correlation (7.26) and 
/3 = — 1, implies the behavior S3 (r) ~ rlnr at small 



44 



separations r. As seen in figure 43, the graph of S 3 (r) 
in log-log coordinates indeed displays a significant cur- 
vature which is a signature of logarithmic corrections. 
The next step consists in decomposing the velocity in- 
crements 5 r v = v(x + r,t) — v(x,t) into their positive 
5+v and negative 5~v parts. It is clear that 

Sl h %r) = -((5-vf) + ((6+vf), 

S 3 (r) = ((5-vf) + ((5+vf), (7.28) 

so that ({S+v) 3 ) = (5f bs (r) + S , 3 (r))/2. As seen in fig- 
ure 43 the log-log plot of ((<5+u) 3 ) as a function of r 
is nearly a straight line with slope 1.07 very close to 
unity. This observation is confirmed in [91] by indepen- 
dently measuring the PDFs of positive and negative ve- 
locity increments. Assuming that ((S^'v) 3 ') ~ Br, one 
obtains the following prediction for the small-r behav- 
iors of the third-order structure functions 

Sf s (r) ~ -Arlnr + Br, 

S 3 (r] |~ Arlnr + Br. (7.29) 

This suggests that the only difference in the small- 
separation behaviors of S'f bs (r) and Ss(r) is the sign in 
the balance between the leading term a rlnr and the 
subleading term oc r. In a log-log plot this difference 
amounts to shifting the graph away from where it is 
most curved and thus makes it straighter, albeit with a 
(local) slope which is not unity. This explains why signif- 
icant deviations from 1 are observed for (3. Note that a 
similar approach can be used for higher-order structure 
functions. It leads for instance to Si(r) s=s Cr — Dr 4/>3 , 
where C and D are two positive constants. The nega- 
tive sign before the sub- leading term (r 4 / 3 ) is crucial. 
It implies that, for any finite r, a naive power-law fit 
to S4 can yield a scaling exponent less than unity. The 
presence of sub-leading, power-law terms with oppo- 
site signs also explains the small apparent "anomalous" 
scaling behavior observed for other values of p in the 
simulations. Note that similar artifacts involving two 
competing power-laws have been described in [16,7]. 

The work reported in this section indicates that a naive 
interpretation of numerical measurements might result 
in predicting artificial anomalous scaling laws. In the 
case of Burgers turbulence for which high-resolution nu- 
merics are available and statistical convergence of the 
averages can be guaranteed, we have seen that it is not 
too difficult to identify the numerical artifacts which 
are responsible for such a masquerading. However this 
is not always the case. For instance, it seems reason- 
able enough to claim that attacking the problem of mul- 
tiscaling in spatially extended nonlinear systems, such 
as Navier-Stokes turbulence, requires considerable the- 
oretical insight that must supplement sophisticated and 
heavy numerical simulations and experiments. Note fi- 
nally that, up to now, the question of the presence or not 



of anomalous scaling laws in the Burgers equation with 
a self-similar forcing with exponent — 1 < (3 < remains 
largely open. 

8 Concluding remarks and open questions 

This review summarizes recent work connected with the 
Burgers equation. Originally this model was introduced 
as a simplification of the Navier-Stokes equation with 
the hope of shedding some light on issues such as turbu- 
lence. This hope did not materialize. Nevertheless many 
of the interesting questions that have been addressed 
for Burgers turbulence are eventually transpositions of 
similar questions for Navier-Stokes turbulence. One par- 
ticularly important instance is the issue of universality 
with respect to the form of the forcing and of the initial 
condition. For Burgers turbulence most of the universal 
features, such as scaling exponents or functional forms 
of PDF tails are dominated by the presence of shocks 
and other singularities in the solution. This applies both 
to the case of decaying turbulence driven by random ini- 
tial conditions and randomly forced turbulence. In the 
latter case one is mostly interested in analysis of station- 
ary properties of solutions, for example stationary dis- 
tribution for velocity increments or gradients. Another 
set of questions is motivated by more mathematical con- 
siderations. It mainly concerns the construction of a sta- 
tionary invariant measure when Burgers dynamics in a 
finite-size domain is supplemented by an external ran- 
dom source of energy. Again it has been shown that the 
presence of shocks, and in particular of global shocks, 
plays a crucial role in the construction of the statisti- 
cally stationary solution. Both physical and mathemat- 
ical questions lead to a similar answer: one first needs to 
describe and control shocks. The main message to retain 
for hydrodynamical turbulence is hence a strong confir- 
mation of the common wisdom that it cannot be fully 
understood without a detailed description of singulari- 
ties. Moreover, the behavior depends not only on the lo- 
cal structure of singularities, but also on their distribu- 
tion at larger scales. Here a word of caution: for incom- 
pressible fully developped Navier-Stokes turbulence, we 
have no evidence that the universal scaling properties 
observed in experiments and simulations stem from real 
singularities. Indeed the issue of a finite-time blow-up of 
the three-dimensional Eulcr equation is still open (see, 
e.g. [56]). Another important observation that can be 
drawn from the study of Burgers turbulence is that both 
the tools used and the answers obtained strongly de- 
pend on the kind of setting one considers: decay versus 
forced turbulence, finite-size versus infinite-size systems, 
smooth versus self-similar forcing, etc. 

Besides turbulence, the random Burgers equation has 
various applications in cosmology, in non-equilibrium 
statistical physics and in disordered media. Among 
them, the connection to the problem of directed poly- 
mers has attracted much attention. As already noted in 
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the Introduction, there is a mathematical equivalence 
between the zero-viscosity limit of the forced Burgers 
equation and the zero-temperature limit for directed 
polymers. We have seen in section 5.4 that the so-called 
KPZ scaling, which usually is derived for a finite tem- 
perature, can be established can be established also in 
the zero-temperatur limit, using the action minimizcr 
representation. Such an observation leads to two related 
questions: to what extent can the limit of zero temper- 
ature give an insight into finite-temperature polymer 
dynamics and how can the global minimizer formalism 
be extended to tackle the finite-temperature setting? 
It looks plausible that in polymer dynamics, or more 
generally in the study of random walks in a random 
potential, the trajectories carrying most of the Gibbs 
probability weight are defining corridors in space time. 
These objects can concentrate near the trajectories of 
global minimizers but, at the moment, there is no for- 
malism to describe them, nor attempts to quantify their 
contribution to the Gibbs statistics. 

Another important open question concerns the multi- 
dimensional extensions of the Burgers equation. As we 
have seen, when the forcing is potential, the potential 
character of the velocity field is conserved by the dynam- 
ics. This leads to the construction of stationary solutions 
which carry many similarities with the one-dimensional 
case. Up to now there is only limited understanding of 
what happens when the potentiality assumption of the 
flow is dropped. This problem has of course concrete 
applications in gas dynamics and for disperse inelastic 
granular media (see, e.g., [12]). An interesting question 
concerns the construction of the limit of vanishing vis- 
cosity, given that the Hopf-Cole transformation is in- 
applicable in the non-potential case. Understanding ex- 
tensions of the viscous limiting procedure to the non- 
potential case might give new insight into the problem of 
the large Reynolds number limit in incompressible tur- 
bulence. Another question related to non-potential flow 
concerns the interactions between vorticity and shocks. 
For instance, in two dimensions the vorticity is trans- 
ported by the flow. This results in its growth in the highly 
compressible regions of the flow. The various singulari- 
ties of the velocity field should hence be strongly affected 
by the flow rotation and, in particular, the shocks are 
expected to have a spiraling structure. 

We finish with few remarks on open mathematical prob- 
lems. As we have seen in the one-dimensional case, one 
can rigorously prove hyperbolicity of the global mini- 
mizer. In the multi-dimensional case it is also possible 
to establish the existence and, in many cases, unique- 
ness of the global minimizer. However, the very impor- 
tant question of its hyperbolicity is still an open prob- 
lem. If proven, hyperbolicity would allow for rigorous 
analysis of the regularity properties of the stationary so- 
lutions and of the topological shocks. There are many 
interesting problems - even basic issues of existence and 
uniqueness - in the non-compact case where at present a 



mathematical theory is basically absent. Finally, a very 
challenging open problem concerns the extension of the 
results on the evolution of matter inside shocks to the 
case of general Hamilton- Jacobi equations. 
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